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

Add support for number of connections different from number of faces.

Browse Source 
This is a precursor to supporting wells as faces.
This commit is contained in:
Bård Skaflestad
2010-08-09 08:51:01 +00:00
parent cd829252ca
commit b196437776
4 changed files with 175 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
mex_ip_simple.c
View File

@@ -1,4 +1,6 @@
#include <string.h>
#include <mex.h>
#include "mrst_api.h"
#include "mimetic.h"
@@ -76,23 +78,30 @@ verify_structural_consistency(int nlhs, int nrhs, const mxArray *prhs[])
{
int rock_ok;
int grid_ok;
int conn_ok;
mxAssert((nlhs == 1) && (nrhs == 2),
"Must be called with exactly two inputs and one output.");
mxAssert((nlhs == 1) && (nrhs == 4),
"Must be called with exactly four inputs and one output.");
mxAssert(mxIsStruct(prhs[0]) && mxIsStruct(prhs[1]),
"Inputs must be STRUCTs.");
"First two inputs must be STRUCTs.");
mxAssert((mxGetNumberOfElements(prhs[0]) == 1) &&
(mxGetNumberOfElements(prhs[1]) == 1),
"Inputs must be 1-by-1 STRUCT arrays.");
"First two inputs must be 1-by-1 STRUCT arrays.");
grid_ok = verify_grid_structure(prhs[0]);
/* rock must contain a field 'perm' */
rock_ok = mxGetFieldNumber(prhs[1], "perm") >= 0;
return grid_ok && rock_ok;
/* connection structure */
conn_ok = (mxIsDouble(prhs[2]) || mxIsInt32(prhs[2])) &&
(mxIsDouble(prhs[3]) || mxIsInt32(prhs[3])) &&
(mxGetNumberOfElements(prhs[2]) <
mxGetNumberOfElements(prhs[3]));
return grid_ok && rock_ok && conn_ok;
}
@@ -122,11 +131,72 @@ deallocate_face_data(int *fneighbour, double *fnormal, double *fcentroid)
}
/* ------------------------------------------------------------------ */
static void
copy_int_vector(const mxArray *M_v, int *v)
/* ------------------------------------------------------------------ */
{
size_t n, i;
int *pi;
double *pd;
n = mxGetNumberOfElements(M_v);
if (mxIsDouble(M_v)) {
pd = mxGetPr(M_v);
for (i = 0; i < n; i++) {
mxAssert ((INT_MIN <= pd[i]) && (pd[i] <= INT_MAX),
"Data element outside range for INT32");
v[i] = pd[i];
}
} else {
pi = mxGetData(M_v);
memcpy(v, pi, n * sizeof *v);
}
}
/* ------------------------------------------------------------------ */
static void
extract_connection_data(const mxArray *M_nconn, const mxArray *M_conn,
const int nc, int **nconn, int **conn,
int *sum_nconn, int *sum_nconn2)
/* ------------------------------------------------------------------ */
{
size_t i, n;
n = mxGetNumberOfElements(M_conn);
*nconn = mxMalloc(nc * sizeof *nconn);
*conn = mxMalloc(n * sizeof *conn );
copy_int_vector(M_nconn, *nconn);
copy_int_vector(M_conn, *conn );
/* Adjust connections for 1-based indexing. */
for (i = 0; i < n; i++) {
(*conn)[i] -= 1;
}
/* Count connections */
n = nc;
*sum_nconn = *sum_nconn2 = 0;
for (i = 0; i < n; i++) {
*sum_nconn += (*nconn)[i];
*sum_nconn2 += (*nconn)[i] * (*nconn)[i];
}
}
/* ------------------------------------------------------------------ */
static int
extract_cell_data(const mxArray *G, int d,
int *max_ncf, int *sum_ncf, int *sum_ncf2,
int **ncfaces, int **cfaces,
extract_cell_data(const mxArray *G,
const mxArray *M_nconn, const mxArray *M_conn,
int d, int *max_ncf, int *sum_nconn, int *sum_nconn2,
int **ncfaces, int **nconn, int **conn,
double **ccentroids, double **cvolumes)
/* ------------------------------------------------------------------ */
{
@@ -134,7 +204,10 @@ extract_cell_data(const mxArray *G, int d,
int ncells = getNumberOfCells(G);
int i;
*max_ncf = *sum_ncf = *sum_ncf2 = 0;
extract_connection_data(M_nconn, M_conn, ncells,
nconn, conn, sum_nconn, sum_nconn2);
*max_ncf = 0;
for(i=0; i<ncells; ++i)
{
@@ -142,11 +215,8 @@ extract_cell_data(const mxArray *G, int d,
(*ncfaces)[i] = n;
*max_ncf = MAX(*max_ncf, n);
*sum_ncf += n;
*sum_ncf2 += n*n;
}
*cfaces = getCellFaces(G);
*ccentroids = getCellCentroids(G);
*cvolumes = getCellVolumes(G);
@@ -156,11 +226,12 @@ extract_cell_data(const mxArray *G, int d,
/* ------------------------------------------------------------------ */
static void
deallocate_cell_data(int *ncfaces, int *cfaces, double *ccentroids)
deallocate_cell_data(int *ncfaces, int *nconn, int *conn, double *ccentroids)
/* ------------------------------------------------------------------ */
{
mxFree(ccentroids);
mxFree(cfaces);
mxFree(conn);
mxFree(nconn);
mxFree(ncfaces);
}
@@ -175,7 +246,7 @@ deallocate_perm_data(double *K)
/*
* BI = mex_ip_simple(G, rock)
* BI = mex_ip_simple(G, rock, nconn, conn)
*/
/* ------------------------------------------------------------------ */
@@ -189,7 +260,7 @@ mexFunction(int nlhs, mxArray *plhs[],
int *fneighbour;
double *farea, *fnormal, *fcentroid;
int *ncfaces, *cfaces;
int *ncfaces, *nconn, *conn;
double *ccentroids, *cvolumes;
double *perm;
@@ -203,25 +274,26 @@ mexFunction(int nlhs, mxArray *plhs[],
nfaces = extract_face_data(prhs[0], d,
&fneighbour, &farea, &fnormal, &fcentroid);
ncells = extract_cell_data(prhs[0], d,
ncells = extract_cell_data(prhs[0], prhs[2], prhs[3], d,
&max_ncf, &sum_ncf, &sum_ncf2, &ncfaces,
&cfaces, &ccentroids, &cvolumes);
&nconn, &conn, &ccentroids, &cvolumes);
perm = getPermeability(mxGetField(prhs[1], 0, "perm"), d);
if ((ncells > 0) && (nfaces > 0)) {
plhs[0] = mxCreateDoubleMatrix(sum_ncf2, 1, mxREAL);
Binv = mxGetPr(plhs[0]);
Binv = mxGetPr(plhs[0]); /* mxCreateDoubleMatrix == ZEROS */
mim_ip_simple_all(ncells, d, max_ncf, ncfaces, cfaces, fneighbour,
fcentroid, fnormal, farea, ccentroids, cvolumes,
perm, Binv);
/* mim_ip_simple_all() requires zeroed target array... */
mim_ip_simple_all(ncells, d, max_ncf, ncfaces, nconn, conn,
fneighbour, fcentroid, fnormal, farea,
ccentroids, cvolumes, perm, Binv);
} else {
plhs[0] = mxCreateDoubleScalar(mxGetNaN());
}
deallocate_perm_data(perm);
deallocate_cell_data(ncfaces, cfaces, ccentroids);
deallocate_cell_data(ncfaces, nconn, conn, ccentroids);
deallocate_face_data(fneighbour, fnormal, fcentroid);
} else {
plhs[0] = mxCreateDoubleScalar(mxGetNaN());

26
mex_ip_simple.m
View File

@@ -2,16 +2,24 @@ function varargout = mex_ip_simple(varargin)
%Compute 'ip_simple' inner product values using compiled C code.
%
% SYNOPSIS:
% BI = mex_ip_simple(G, rock)
% BI = mex_ip_simple(G, rock, nconn, conn)
%
% 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.
%
% RETURNS:
% BI - A SUM(DIFF(G.cells.facePos) .^ 2)-by-1 array of inner product
% values, ordered by the cells of the input grid.
% BI - A SUM(nconn .^ 2)-by-1 array of inner product values, ordered by
% the cells of the input grid.
%
% NOTE:
% As the return value 'BI' is but a simple data array value, it must be
@@ -30,12 +38,14 @@ 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);
% BI = mex_ip_simple(G, rock, nconn, conn);
% toc(t0)
%
% [i, j] = blockDiagIndex(diff(G.cells.facePos), ...
% diff(G.cells.facePos));
% [i, j] = blockDiagIndex(nconn, nconn);
%
% S = struct('BI', sparse(i, j, BI), 'type', 'hybrid', 'ip', 'ip_simple')
%

76
mimetic.c
View File

@@ -1,9 +1,9 @@
#include <assert.h>
#include <stddef.h>
#include <stdlib.h>
#if defined(COMPILING_FOR_MATLAB) && COMPILING_FOR_MATLAB
#include "mex.h"
#include "matrix.h"
#include <mex.h>
#define MAT_SIZE_T mwSignedIndex
#endif
@@ -16,8 +16,9 @@
/* ------------------------------------------------------------------ */
void
mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
int *fneighbour, double *fcentroid, double *fnormal,
mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
int *nconn, int *cf,
int *fneighbour, double *fcentroid, double *fnormal,
double *farea, double *ccentroid, double *cvol,
double *perm, double *Binv)
/* ------------------------------------------------------------------ */
@@ -56,11 +57,11 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
}
}
mim_ip_simple(nf, d, cvol[c], A, &perm[c * d * d],
C, N, &Binv[fpos2], work, lwork);
mim_ip_simple(nf, nconn[c], d, cvol[c], &perm[c * d * d],
C, A, N, &Binv[fpos2], work, lwork);
fpos += nf;
fpos2 += nf * nf;
fpos += nconn[c];
fpos2 += nconn[c] * nconn[c];
}
}
@@ -70,32 +71,41 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
/* ------------------------------------------------------------------ */
void
mim_ip_simple(int nf, int d, double vol,
double *A, double *K, double *C, double *N,
double *Binv, double *work, int lwork)
mim_ip_simple(int nf, int nconn, int d,
double v, double *K, double *C,
double *A, double *N,
double *Binv,
double *work, int lwork)
/* ------------------------------------------------------------------ */
{
mim_ip_span_nullspace(nf, d, C, A, Binv, work, lwork);
mim_ip_linpress_exact(nf, d, vol, N, K, Binv, C);
mim_ip_span_nullspace(nf, nconn, d, C, A, Binv, work, lwork);
mim_ip_linpress_exact(nf, nconn, d, v, K, N, Binv, work, lwork);
}
/* ------------------------------------------------------------------ */
void
mim_ip_span_nullspace(int nf, int d, double *C, double *A,
double *X, double *work, int nwork)
mim_ip_span_nullspace(int nf, int nconn, int d,
double *C,
double *A,
double *X,
double *work, int nwork)
/* ------------------------------------------------------------------ */
{
MAT_SIZE_T m, n, k, ld, info, lwork;
MAT_SIZE_T m, n, k, ldC, ldX, info, lwork;
int i, j;
double a1, a2;
double tau[3] = { 0.0 }; /* No more than 3 spatial dimensions */
/* Step 1) X <- I */
for (i = 0; i < nf*nf; i++) { X[i ] = 0.0; }
for (i = 0; i < nf ; i++) { X[i * (nf + 1)] = 1.0; }
/* Step 1) X(1:nf, 1:nf) <- I_{nf} */
for (j = 0; j < nf; j++) {
for (i = 0; i < nf; i++) {
X[i + j*nconn] = 0.0;
}
X[j * (nconn + 1)] = 1.0;
}
/* Step 2) C <- orth(A * C) */
for (j = 0; j < d; j++) {
@@ -104,24 +114,25 @@ mim_ip_span_nullspace(int nf, int d, double *C, double *A,
}
}
m = nf; n = d; ld = nf; k = d; lwork = nwork;
dgeqrf_(&m, &n, C, &ld, tau, work, &lwork, &info);
dorgqr_(&m, &n, &k, C, &ld, tau, work, &lwork, &info);
m = nf; n = d; ldC = nf; k = d; lwork = nwork;
dgeqrf_(&m, &n, C, &ldC, tau, work, &lwork, &info);
dorgqr_(&m, &n, &k, C, &ldC, tau, work, &lwork, &info);
/* Step 3) X <- A * (X - C*C') * A */
ldX = nconn;
a1 = -1.0; a2 = 1.0;
dsyrk_("Upper Triangular", "No Transpose",
&m, &n, &a1, C, &ld, &a2, X, &ld);
&m, &n, &a1, C, &ldC, &a2, X, &ldX);
for (j = 0; j < nf; j++) {
for (i = 0; i <= j; i++) {
X[i + j*nf] *= A[i] * A[j];
X[i + j*nconn] *= A[i] * A[j];
}
}
/* Account for DSYRK only assigning upper triangular part. */
for (j = 0; j < nf; j++) {
for (i = j + 1; i < nf; i++) {
X[i + j*nf] = X[j + i*nf];
X[i + j*nconn] = X[j + i*nconn];
}
}
}
@@ -129,14 +140,19 @@ mim_ip_span_nullspace(int nf, int d, double *C, double *A,
/* ------------------------------------------------------------------ */
void
mim_ip_linpress_exact(int nf, int d, double vol, double *N, double *K,
double *Binv, double *T)
mim_ip_linpress_exact(int nf, int nconn, int d,
double vol, double *K,
double *N,
double *Binv,
double *work, int lwork)
/* ------------------------------------------------------------------ */
{
MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
int i;
double a1, a2, t;
assert (lwork >= d * nf);
t = 0.0;
for (i = 0; i < d; i++) {
t += K[i + i*d];
@@ -147,12 +163,12 @@ mim_ip_linpress_exact(int nf, int d, double vol, double *N, double *K,
ld1 = nf ; ld2 = d ;
a1 = 1.0; a2 = 0.0;
dgemm_("No Transpose", "No Transpose", &m, &n, &k,
&a1, N, &ld1, K, &ld2, &a2, T, &ld1);
&a1, N, &ld1, K, &ld2, &a2, work, &ld1);
/* Step 5) Binv <- (N*K*N' + t*X) / vol */
a1 = 1.0 / vol ;
a2 = 6.0 * t / (d * vol);
ldBinv = nf;
ldBinv = nconn;
dgemm_("No Transpose", "Transpose", &m, &m, &n,
&a1, T, &ld1, N, &ld1, &a2, Binv, &ldBinv);
&a1, work, &ld1, N, &ld1, &a2, Binv, &ldBinv);
}

23
mimetic.h
View File

@@ -2,17 +2,26 @@
#define MIMETIC_H_INCLUDED
void mim_ip_span_nullspace(int nf, int d, double *C, double *A,
double *X, double *work, int lwork);
void mim_ip_span_nullspace(int nf, int nconn, int d,
double *C,
double *A,
double *X,
double *work, int lwork);
void mim_ip_linpress_exact(int nf, int d, double vol, double *N,
double *K, double *Binv, double *T);
void mim_ip_linpress_exact(int nf, int nconn, int d,
double vol, double *K,
double *N,
double *Binv,
double *work, int lwork);
void mim_ip_simple(int nf, int d, double v, double *A, double *K,
double *C, double *N, double *Binv,
void mim_ip_simple(int nf, int nconn, int d,
double v, double *K, double *C,
double *A, double *N,
double *Binv,
double *work, int lwork);
void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
int *nconn, int *conn,
int *fneighbour, double *fcentroid, double *fnormal,
double *farea, double *ccentroid, double *cvol,
double *perm, double *Binv);