diff --git a/mex_compute_coarse_contrib.c b/mex_compute_coarse_contrib.c new file mode 100644 index 00000000..b004530a --- /dev/null +++ b/mex_compute_coarse_contrib.c @@ -0,0 +1,230 @@ +#include +#include + +#include + +#include "partition.h" +#include "coarse_sys.h" + + +/* ---------------------------------------------------------------------- */ +static int +args_ok(int nlhs, int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int i, ok; + + ok = (nlhs == 2) && (nrhs == 6); + ok = ok && mxIsDouble(prhs[0]); + ok = ok && mxIsDouble(prhs[1]); + + for (i = 2; ok && (i < nrhs); i++) { + ok = mxIsDouble(prhs[i]) || mxIsInt32(prhs[i]); + } + + return ok; +} + + +/* ---------------------------------------------------------------------- */ +static void +extract_int_vec(const mxArray *M_v, int *v) +/* ---------------------------------------------------------------------- */ +{ + size_t i, n; + + int *pi; + double *pd; + + n = mxGetNumberOfElements(M_v); + + if (mxIsDouble(M_v)) { + pd = mxGetPr(M_v); + + for (i = 0; i < n; i++) { v[i] = pd[i]; } + } else { + pi = mxGetData(M_v); + + for (i = 0; i < n; i++) { v[i] = pi[i]; } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +add_constant(size_t n, int a, int *v) +/* ---------------------------------------------------------------------- */ +{ + size_t i; + + for (i = 0; i < n; i++) { v[i] += a; } +} + + +/* ---------------------------------------------------------------------- */ +static int +max_element(size_t n, int *v) +/* ---------------------------------------------------------------------- */ +{ + int ret; + size_t i; + + ret = v[0]; + for (i = 1; i < n; i++) { + ret = (v[i] > ret) ? v[i] : ret; + } + + return ret; +} + + +/* ---------------------------------------------------------------------- */ +static int +max_diff(size_t n, int *p) +/* ---------------------------------------------------------------------- */ +{ + int d, ret; + size_t i; + + ret = -1; + for (i = 0; i < n; i++) { + d = p[i + 1] - p[i]; + ret = (d > ret) ? d : ret; + } + + return ret; +} + + +/* ---------------------------------------------------------------------- */ +static int +extract_partition(const mxArray *M_p, int *nc, int *nb, int **p) +/* ---------------------------------------------------------------------- */ +{ + int ret; + + *nc = mxGetNumberOfElements(M_p); + *p = mxMalloc(*nc * sizeof **p); + + if (*p != NULL) { + extract_int_vec(M_p, *p); + add_constant(*nc, -1, *p); + + *nb = max_element(*nc, *p) + 1; + + ret = 1; + } else { + ret = 0; + } + + return ret; +} + + +/* + * [cell_ip, Binv] = mex_compute_coarse_contrib(BIf, Psi, p, pconn, ... + * dof_pos, blk_ncf) + */ +/* ---------------------------------------------------------------------- */ +void +mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok, nc, nb, b, i, tmp, max_bcells, max_nconn, max_nbpair; + int *p, *pconn, *blk_ncf, *b2c_pos, *b2c; + + size_t totblkdof; + + char errmsg[1023 + 1]; + + double *BI, *work, *totmob; + + struct coarse_sys sys; + + if (args_ok(nlhs, nrhs, prhs)) { + BI = mxGetPr(prhs[0]); + sys.basis = mxGetPr(prhs[1]); + ok = extract_partition(prhs[2], &nc, &nb, &p); + + if (ok) { + ok = partition_allocate_inverse(nc, nb - 1, + &b2c_pos, &b2c); + + pconn = mxMalloc((nc + 1) * sizeof *pconn ); + blk_ncf = mxMalloc( nb * sizeof *blk_ncf ); + sys.dof_pos = mxMalloc((nb + 1) * sizeof *sys.dof_pos); + sys.bf_pos = mxCalloc( nb + 1 , sizeof *sys.bf_pos ); + sys.ip_pos = mxCalloc( nb + 1 , sizeof *sys.ip_pos ); + + if (ok && (pconn != NULL) && (blk_ncf != NULL) && + (sys.dof_pos != NULL) && (sys.bf_pos != NULL) && + (sys.ip_pos != NULL)) { + partition_invert(nc, p, b2c_pos, b2c); + + extract_int_vec(prhs[3], pconn ); + extract_int_vec(prhs[4], sys.dof_pos); + extract_int_vec(prhs[5], blk_ncf ); + + add_constant(nc + 1, -1, pconn); + add_constant(nb + 1, -1, sys.dof_pos); + + max_nconn = max_diff(nc, pconn); + max_bcells = max_diff(nb, b2c_pos); + max_nbpair = max_diff(nb, sys.dof_pos); + max_nbpair = max_nbpair * (max_nbpair + 1) / 2; + + totblkdof = 0; + for (b = 1; b <= nb; b++) { + tmp = sys.dof_pos[b] - sys.dof_pos[b - 1]; + + totblkdof += tmp * tmp; + + sys.bf_pos[b] = sys.bf_pos[b - 1] + tmp*blk_ncf[b - 1]; + + tmp = tmp * (tmp + 1) / 2; + sys.ip_pos[b] = sys.ip_pos[b - 1] + tmp*(b2c_pos[b] - b2c_pos[b - 1]); + } + + plhs[0] = mxCreateDoubleMatrix(sys.ip_pos[nb], 1, mxREAL); + plhs[1] = mxCreateDoubleMatrix(totblkdof , 1, mxREAL); + work = mxMalloc((max_bcells + max_nbpair*(max_nbpair + 1)/2) * + sizeof *work); + totmob = mxMalloc(nc * sizeof *totmob); + + if ((plhs[0] != NULL) && (plhs[1] != NULL) && + (work != NULL) && (totmob != NULL)) { + sys.cell_ip = mxGetPr(plhs[0]); + sys.Binv = mxGetPr(plhs[1]); + + for (i = 0; i < nc; i++) { totmob[i] = 1.0; } + + coarse_sys_compute_cell_ip(nc, max_nconn, nb, pconn, BI, + b2c_pos, b2c, &sys); + coarse_sys_compute_Binv(nb, max_bcells, totmob, + b2c_pos, b2c, &sys, work); + } + + if (totmob != NULL) { mxFree(totmob); } + if (work != NULL) { mxFree(work); } + } + + if (sys.ip_pos != NULL) { mxFree(sys.ip_pos); } + if (sys.bf_pos != NULL) { mxFree(sys.bf_pos); } + if (sys.dof_pos != NULL) { mxFree(sys.dof_pos); } + if (blk_ncf != NULL) { mxFree(blk_ncf); } + if (pconn != NULL) { mxFree(pconn); } + + partition_deallocate_inverse(b2c_pos, b2c); + } + + if (p != NULL) { mxFree(p); } + } else { + sprintf(errmsg, + "Calling sequence is\n\t" + "[cell_ip, Binv] = %s(BIf, Psi, p, pconn, dof_pos, blk_ncf)", + mexFunctionName()); + + mexErrMsgTxt(errmsg); + } +} diff --git a/test_cs_assembly.m b/test_cs_assembly.m new file mode 100644 index 00000000..47e140ae --- /dev/null +++ b/test_cs_assembly.m @@ -0,0 +1,57 @@ +g = computeGeometry(cartGrid([4, 4])); +rock.perm = ones([g.cells.num, 1]); + +[BI, pconn, conn] = mex_ip_simple(g, rock); + +nconn = double(diff(pconn)); +[i, j] = blockDiagIndex(nconn, nconn); + +S = struct('BI', sparse(i, j, BI), 'type', 'hybrid', 'ip', 'ip_simple'); +x = initResSol(g, 0); + + +fluid = initSingleFluid('mu', 1, 'rho', 1); +mu = fluid.properties(x); +s = fluid.saturation(x); +kr = fluid.relperm(s, x); + +mob = sum(bsxfun(@rdivide, kr, mu), 2); + +p = mex_partition_ui(double(g.cells.indexMap), g.cartDims, [2, 2]); + +[b2c_pos, b2c] = mex_partition_invert(p); + +cg = mex_generate_coarsegrid(g, p, 5); +cs = generateCoarseSystem(g, rock, S, cg, mob); + +BPsi = basisMatrixHybrid(g, cg, cs); +Psi = S.BI * BPsi; + + +cellno = rldecode(1:g.cells.num, diff(g.cells.facePos), 2) .'; +a = sortrows([double(p(cellno)), (1:size(g.cells.faces,1))']); +blkcf_pos = cumsum([1; accumarray(a(:,1), 1)]); + + +blkno = rldecode(1:cg.cells.num, diff(cg.cells.facePos), 2) .'; +t = [blkno(cs.activeCellFaces), (1 : numel(cs.activeCellFaces)) .']; +tpos = cumsum([1; accumarray(t(:,1), 1)]); + +Psi_vals = []; +for b = 1 : numel(tpos) - 1, + for cf = reshape(t(tpos(b) : tpos(b+1) - 1, 2), 1, []), + Psi_vals = [Psi_vals; ... + full(Psi(a(blkcf_pos(b) : blkcf_pos(b+1) - 1, 2), cf))]; + end +end + + +dof = zeros([cg.faces.num, 1]); +dof(cs.activeFaces) = 1 : numel(cs.activeFaces); + + +[cell_ip, Binv] = mex_compute_coarse_contrib(BI, Psi_vals, p, pconn, ... + tpos, diff(blkcf_pos)); + +[ss, rr, ff, ll] = ... + mex_schur_comp_symm(Binv, tpos, dof(cg.cells.faces(cs.activeCellFaces)));