From 704d7eab5c63baa14235b353b4bdb6000ebb6e1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:13:41 +0200 Subject: [PATCH] Backed out changeset 5e3d75476d64 In preparation of restoring backwards compatibility in 'cfs_tpfa' module. --- src/cfs_tpfa.c | 841 +++++++++++++++++++++++++++++-------------------- src/cfs_tpfa.h | 18 +- 2 files changed, 517 insertions(+), 342 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 8ce3efc0b..0fea19f90 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -24,25 +24,33 @@ struct densrat_util { MAT_SIZE_T *ipiv; - - double residual; double *lu; - double *t1; - double *t2; - double *mat_row; - double *coeff; - double *linsolve_buffer; + + double *x; + double *Ai_y; + double *psum; + + /* Storage */ + double *ddata; }; struct cfs_tpfa_impl { - int is_incomp; + /* Reservoir flow */ + double *ctrans; + double *accum; /* One entry per component per face */ - double *compflux_f; /* A_{ij} v_{ij} */ - double *compflux_deriv_f; /* A_{ij} \partial_{p} v_{ij} */ + double *masstrans_f; /* RB^{-1} [ phase-mobility ] */ + double *gravtrans_f; /* RB^{-1} [ grav + capillary ] */ - double *flux_work; + /* Compressible completion flow definition */ + double *wtrans; /* WI * sum((Ac \ Ap) [ pmob ]) */ + double *wgpot; /* WI * sum((Ac \ Ap) [ pmob ]'*WdP) */ + + /* One entry per component per completion/perforation */ + double *masstrans_p; /* RB^{-1} [ phase-mobility ] */ + double *gravtrans_p; /* RB^{-1} [ grav ] */ /* Scratch array for face pressure calculation */ double *scratch_f; @@ -60,7 +68,7 @@ deallocate_densrat(struct densrat_util *ratio) /* ---------------------------------------------------------------------- */ { if (ratio != NULL) { - free(ratio->lu); + free(ratio->ddata); free(ratio->ipiv); } @@ -70,41 +78,46 @@ deallocate_densrat(struct densrat_util *ratio) /* ---------------------------------------------------------------------- */ static struct densrat_util * -allocate_densrat(size_t max_conn, int np) +allocate_densrat(grid_t *g, well_t *w, int np) /* ---------------------------------------------------------------------- */ { - size_t alloc_sz, n_buffer_col; - struct densrat_util *ratio; + int ntotperf; + size_t nglobconn, ntotconn, ddata_sz; - ratio = malloc(1 * sizeof *ratio); + struct densrat_util *new; - if (ratio != NULL) { - n_buffer_col = 1; /* z */ - n_buffer_col += 1 * max_conn; /* A_{ij} v_{ij} */ - n_buffer_col += 2 * max_conn; /* A_{ij} \partial_{p} v_{ij} */ + new = malloc(1 * sizeof *new); - alloc_sz = np * np; /* lu */ - alloc_sz += 2 * np; /* t1, t2 */ - alloc_sz += max_conn * 1 ; /* mat_row */ - alloc_sz += (max_conn + 1) * 1 ; /* coeff */ - alloc_sz += n_buffer_col * np; /* linsolve_buffer */ - - ratio->ipiv = malloc(np * sizeof *ratio->ipiv); - ratio->lu = malloc(alloc_sz * sizeof *ratio->lu ); - - if ((ratio->ipiv == NULL) || (ratio->lu == NULL)) { - deallocate_densrat(ratio); - ratio = NULL; + if (new != NULL) { + if (w != NULL) { + ntotperf = w->well_connpos[ w->number_of_wells ]; } else { - ratio->t1 = ratio->lu + (np * np); - ratio->t2 = ratio->t1 + (1 * np); - ratio->mat_row = ratio->t2 + (1 * np); - ratio->coeff = ratio->mat_row + (max_conn * 1 ); - ratio->linsolve_buffer = ratio->coeff + ((max_conn + 1) * 1 ); + ntotperf = 0; + } + + nglobconn = MAX(g->number_of_faces , ntotperf); + ntotconn = MAX(g->cell_facepos[ g->number_of_cells ], ntotperf); + + ddata_sz = np * np; /* lu */ + ddata_sz += np * nglobconn; /* x */ + ddata_sz += np * ntotconn; /* Ai_y */ + ddata_sz += 1 * ntotconn; /* psum */ + + new->ipiv = malloc(np * sizeof *new->ipiv); + new->ddata = malloc(ddata_sz * sizeof *new->ddata); + + if ((new->ipiv == NULL) || (new->ddata == NULL)) { + deallocate_densrat(new); + new = NULL; + } else { + new->lu = new->ddata + 0 ; + new->x = new->lu + np * np ; + new->Ai_y = new->x + np * nglobconn; + new->psum = new->Ai_y + np * ntotconn ; } } - return ratio; + return new; } @@ -124,7 +137,7 @@ impl_deallocate(struct cfs_tpfa_impl *pimpl) /* ---------------------------------------------------------------------- */ static struct cfs_tpfa_impl * -impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) +impl_allocate(grid_t *G, well_t *W, int np) /* ---------------------------------------------------------------------- */ { size_t nnu, ngconn, nwperf; @@ -145,18 +158,24 @@ impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) ddata_sz = 2 * nnu; /* b, x */ /* Reservoir */ - ddata_sz += np * G->number_of_faces ; /* compflux_f */ - ddata_sz += np * (2 * G->number_of_faces); /* compflux_deriv_f */ + ddata_sz += 1 * ngconn; /* ctrans */ + ddata_sz += 1 * G->number_of_cells; /* accum */ + ddata_sz += np * G->number_of_faces; /* masstrans_f */ + ddata_sz += np * G->number_of_faces; /* gravtrans_f */ - ddata_sz += np * (1 + 2) ; /* flux_work */ + /* Wells */ + ddata_sz += 1 * nwperf; /* wtrans */ + ddata_sz += 1 * nwperf; /* wgpot */ + ddata_sz += np * nwperf; /* masstrans_p */ + ddata_sz += np * nwperf; /* gravtrans_p */ - ddata_sz += 1 * G->number_of_faces ; /* scratch_f */ + ddata_sz += 1 * G->number_of_faces; /* scratch_f */ new = malloc(1 * sizeof *new); if (new != NULL) { new->ddata = malloc(ddata_sz * sizeof *new->ddata); - new->ratio = allocate_densrat(max_conn, np); + new->ratio = allocate_densrat(G, W, np); if (new->ddata == NULL || new->ratio == NULL) { impl_deallocate(new); @@ -260,277 +279,348 @@ construct_matrix(grid_t *G, well_t *W) } +/* ---------------------------------------------------------------------- */ static void -factorise_fluid_matrix(int np, const double *A, struct densrat_util *ratio) +solve_cellsys_core(grid_t *G , + size_t sz , + const double *Ac , + const double *bf , + double *xcf , + double *luAc, + MAT_SIZE_T *ipiv) +/* ---------------------------------------------------------------------- */ { - int np2; - MAT_SIZE_T m, n, ld, info; + int c, i, f; + size_t j, p2; + double *v; - m = n = ld = np; - np2 = np * np; + MAT_SIZE_T nrows, ncols, ldA, ldX, nrhs, info; - memcpy (ratio->lu, A, np2 * sizeof *ratio->lu); - dgetrf_(&m, &n, ratio->lu, &ld, ratio->ipiv, &info); + nrows = ncols = ldA = ldX = sz; + info = 0; - assert (info == 0); + v = xcf; + + for (c = 0, p2 = 0; c < G->number_of_cells; c++) { + /* Define right-hand sides for local systems */ + for (i = G->cell_facepos[c + 0], nrhs = 0; + i < G->cell_facepos[c + 1]; i++, nrhs++) { + f = G->cell_faces[i]; + + for (j = 0; j < sz; j++) { + v[nrhs*sz + j] = bf[f*sz + j]; + } + } + + /* Factor Ac */ + memcpy(luAc, Ac + p2, sz * sz * sizeof *luAc); + dgetrf_(&nrows, &ncols, luAc, &ldA, ipiv, &info); + + /* Solve local systems */ + dgetrs_("No Transpose", &nrows, &nrhs, + luAc, &ldA, ipiv, v, &ldX, &info); + + v += nrhs * sz; + p2 += sz * sz; + } } +/* ---------------------------------------------------------------------- */ static void -solve_linear_systems(int np , - MAT_SIZE_T nrhs , - struct densrat_util *ratio, - double *b ) +small_matvec(size_t n, + int sz, + const double *A, + const double *X, + double *Y) +/* ---------------------------------------------------------------------- */ { - MAT_SIZE_T n, ldA, ldB, info; + size_t i, p1, p2; - n = ldA = ldB = np; + MAT_SIZE_T nrows, ncols, ld, incx, incy; + double a1, a2; - dgetrs_("No Transpose", &n, - &nrhs, ratio->lu, &ldA, ratio->ipiv, - b , &ldB, &info); + nrows = ncols = ld = sz; + incx = incy = 1; - assert (info == 0); -} - - -static void -matvec(int nrow, int ncol, const double *A, const double *x, double *y) -{ - MAT_SIZE_T m, n, ld, incx, incy; - double a1, a2; - - m = ld = nrow; - n = ncol; - incx = incy = 1; - a1 = 1.0; - a2 = 0.0; - - dgemv_("No Transpose", &m, &n, - &a1, A, &ld, x, &incx, - &a2, y, &incy); -} - - -static void -matmat(int np, int ncol, const double *A, const double *B, double *C) -{ - MAT_SIZE_T m, n, k, ldA, ldB, ldC; - double a1, a2; - - m = k = ldA = ldB = ldC = np; - n = ncol; a1 = 1.0; a2 = 0.0; - dgemm_("No Transpose", "No Transpose", &m, &n, &k, - &a1, A, &ldA, B, &ldB, &a2, C, &ldC); + for (i = p1 = p2 = 0; i < n; i++) { + dgemv_("No Transpose", &nrows, &ncols, + &a1, A + p2, &ld, X + p1, &incx, + &a2, Y + p1, &incy); + + p1 += sz; + p2 += sz * sz; + } } +/* ---------------------------------------------------------------------- */ static void -compute_darcyflux_and_deriv(int np, - double trans, - double dp, - const double *pmobf, - const double *gcapf, - double *dflux, - double *dflux_deriv) +solve_cellsys(grid_t *G , + size_t sz, + const double *Ac, + const double *bf, + struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ { - int p; - double a; - - for (p = 0; p < np; p++) { - a = trans * pmobf[p]; - - dflux [ p] = a * (dp + gcapf[p]); - dflux_deriv[0*np + p] = a; /* ignore gravity... */ - dflux_deriv[1*np + p] = - a; - } + solve_cellsys_core(G, sz, Ac, bf, ratio->Ai_y, + ratio->lu, ratio->ipiv); } +/* ---------------------------------------------------------------------- */ static void -compute_compflux_and_deriv(grid_t *G , - int np , - const double *cpress, - const double *trans , - const double *pmobf , - const double *gcapf , - const double *Af , - struct cfs_tpfa_impl *pimpl ) +set_dynamic_trans(grid_t *G , + const double *trans, + struct compr_quantities *cq , + struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ { - int c1, c2, f, np2; - double dp; - double *cflux, *dcflux; + int f, p, i; - np2 = np * np; - - cflux = pimpl->compflux_f; - dcflux = pimpl->compflux_deriv_f; - - for (f = 0; f < G->number_of_faces; - f++, pmobf += np, gcapf += np, Af += np2, - cflux += np, dcflux += 2 * np) { - - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; - - if ((c1 >= 0) && (c2 >= 0)) { - dp = cpress[c1] - cpress[c2]; - - compute_darcyflux_and_deriv(np, trans[f], dp, pmobf, gcapf, - pimpl->flux_work, - pimpl->flux_work + np); - - /* Component flux = Af * v*/ - matvec(np, np, Af, pimpl->flux_work , cflux ); - - /* Derivative = Af * (dv/dp) */ - matmat(np, 2 , Af, pimpl->flux_work + np, dcflux); - } - - /* Boundary connections excluded */ - } -} - - -static int -count_internal_conn(grid_t *G, int c) -{ - int c1, c2, f, i, nconn; - - for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; - - nconn += (c1 >= 0) && (c2 >= 0); - } - - return nconn; -} - - -static int -init_cell_contrib(grid_t *G , - int c , - int np , - double pvol , - double dt , - const double *z , - struct cfs_tpfa_impl *pimpl) -{ - int c1, c2, f, i, conn, nconn; - double *cflx, *dcflx; - - nconn = count_internal_conn(G, c); - - memcpy(pimpl->ratio->linsolve_buffer, z, np * sizeof *z); - - pimpl->ratio->coeff[0] = -pvol; - conn = 1; - - cflx = pimpl->ratio->linsolve_buffer + (1 * np); - dcflx = cflx + (nconn * np); - - for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; c++) { - f = G->cell_faces[i]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; - - if ((c1 >= 0) && (c2 >= 0)) { - memcpy(cflx, pimpl->compflux_f + (f*np + 0), - np * sizeof *cflx); - - memcpy(dcflx, pimpl->compflux_deriv_f + (f*(2 * np) + 0), - 2 * np * sizeof *dcflx); - - cflx += 1 * np; - dcflx += 2 * np; - - pimpl->ratio->coeff[ conn++ ] = dt; + for (f = i = 0; f < G->number_of_faces; f++) { + for (p = 0; p < cq->nphases; p++, i++) { + ratio->x[i] = trans[f] * cq->phasemobf[i]; } } - - assert (conn == nconn + 1); - assert (cflx == pimpl->ratio->linsolve_buffer + (nconn + 1)*np); - - return nconn; } +/* ---------------------------------------------------------------------- */ static void -compute_cell_contrib(grid_t *G , - int c , - int np , - double pvol , - double dt , - const double *z , - const double *Ac , - const double *dAc , - struct cfs_tpfa_impl *pimpl) +set_dynamic_grav(grid_t *G , + flowbc_t *bc , + const double *trans , + const double *gravcap_f, + struct compr_quantities *cq , + struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ { - int c1, c2, f, i, nconn, np2, p; - MAT_SIZE_T nrhs; - double s, dF1, dF2, *dv, *dv1, *dv2; + int f, p, i, c1, c2; + for (f = i = 0; f < G->number_of_faces; f++) { + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; + + if (((c1 >= 0) && (c2 >= 0)) || (bc->type[f] == PRESSURE)) { + for (p = 0; p < cq->nphases; p++, i++) { + ratio->x[i] = trans[f] * gravcap_f[i] * cq->phasemobf[i]; + } + } else { + for (p = 0; p < cq->nphases; p++, i++) { + ratio->x[i] = 0.0; + } + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +set_dynamic_trans_well(well_t *W, + size_t np, + struct completion_data *wdata, + struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ +{ + size_t p, i, k, nconn; + + nconn = W->well_connpos[ W->number_of_wells ]; + + for (i = k = 0; i < nconn; i++) { + for (p = 0; p < np; p++, k++) { + ratio->x[k] = wdata->WI[i] * wdata->phasemob[k]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +set_dynamic_grav_well(well_t *W, + size_t np, + struct completion_data *wdata, + struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ +{ + size_t p, i, k, nconn; + + nconn = W->well_connpos[ W->number_of_wells ]; + + for (i = k = 0; i < nconn; i++) { + for (p = 0; p < np; p++, k++) { + ratio->x[k] = wdata->WI[i] * wdata->gpot[k] * wdata->phasemob[k]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +sum_phase_contrib(grid_t *G , + size_t sz , + const double *xcf, + double *sum) +/* ---------------------------------------------------------------------- */ +{ + int c, i; + size_t j; + + const double *v; + + for (c = i = 0, v = xcf; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++, v += sz) { + + sum[i] = 0.0; + for (j = 0; j < sz; j++) { + sum[i] += v[j]; + } + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_densrat_update(grid_t *G , + struct compr_quantities *cq , + struct densrat_util *ratio, + double *q) +/* ---------------------------------------------------------------------- */ +{ + /* q = Af * x */ + small_matvec(G->number_of_faces, cq->nphases, cq->Af, ratio->x, q); + + /* ratio->Ai_y = Ac \ q */ + solve_cellsys(G, cq->nphases, cq->Ac, q, ratio); + + /* ratio->psum = sum_\alpha ratio->Ai_y */ + sum_phase_contrib(G, cq->nphases, ratio->Ai_y, ratio->psum); +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_densrat_update_well(well_t *W , + struct completion_data *wdata, + struct compr_quantities *cq , + struct densrat_util *ratio, + double *q) +/* ---------------------------------------------------------------------- */ +{ + size_t c, i, nconn, p, np, np2; + MAT_SIZE_T nrows, ncols, ldA, ldX, nrhs, info, incx, incy; + double a1, a2; + + nconn = W->well_connpos[ W->number_of_wells ]; + np = cq->nphases; np2 = np * np; - nconn = init_cell_contrib(G, c, np, pvol, dt, z, pimpl); - nrhs = 1 + (1 + 2)*nconn; /* [z, Af*v, Af*dv] */ + nrows = ncols = ldA = ldX = np; + incx = incy = nrhs = 1; - factorise_fluid_matrix(np, Ac, pimpl->ratio); - solve_linear_systems (np, nrhs, pimpl->ratio, - pimpl->ratio->linsolve_buffer); + a1 = 1.0; + a2 = 0.0; - /* Sum residual contributions over the connections (+ accumulation): - * t1 <- (Ac \ [z, Af*v]) * [-pvol; repmat(dt, [nconn, 1])] */ - matvec(np, nconn + 1, pimpl->ratio->linsolve_buffer, - pimpl->ratio->coeff, pimpl->ratio->t1); + for (i = 0; i < nconn; i++) { + c = W->well_cells[i]; - /* Compute residual in cell 'c' */ - pimpl->ratio->residual = pvol; - for (p = 0; p < np; p++) { - pimpl->ratio->residual += pimpl->ratio->t1[ p ]; - } + /* Compute q = A*x on completion */ + dgemv_("No Transpose", &nrows, &ncols, + &a1, wdata->A + i*np2, &ldA, ratio->x + i*np, &incx, + &a2, q + i*np, &incy); - /* Jacobian row */ - - /* t2 <- A \ ((dA/dp) * t1) */ - matvec(np, np, dAc, pimpl->ratio->t1, pimpl->ratio->t2); - solve_linear_systems(np, 1, pimpl->ratio, pimpl->ratio->t2); - - dF1 = dF2 = 0.0; - for (p = 0; p < np; p++) { - dF1 += pimpl->ratio->t1[ p ]; - dF2 += pimpl->ratio->t2[ p ]; - } - - pimpl->is_incomp = pimpl->is_incomp && (! (fabs(dF2) > 0)); - pimpl->ratio->mat_row[ 0 ] = dF1 - dF2; - - /* Accumulate inter-cell Jacobian contributions */ - dv = pimpl->ratio->linsolve_buffer + (1 + nconn)*np; - for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; i++) { - - f = G->cell_faces[i]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; - - if ((c1 >= 0) && (c2 >= 0)) { - if (c1 == c) { s = 1.0; dv1 = dv + 0 ; dv2 = dv + np; } - else { s = -1.0; dv1 = dv + np; dv2 = dv + 0 ; } - - dF1 = dF2 = 0.0; - for (p = 0; p < np; p++) { - dF1 += dv1[ p ]; - dF2 += dv2[ p ]; - } - - pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; - pimpl->ratio->mat_row[ i ] += s * dt * dF2; + /* Form system RHS */ + for (p = 0; p < np; p++) { + ratio->Ai_y[i*np + p] = q[i*np + p]; } + + /* Factor A in cell 'c' */ + memcpy(ratio->lu, cq->Ac + c*np2, np2 * sizeof *ratio->lu); + dgetrf_(&nrows, &ncols, ratio->lu, &ldA, ratio->ipiv, &info); + + /* Solve local system (=> Ai_y = Ac \ (A*x)) */ + dgetrs_("No Transpose" , &nrows, &nrhs, + ratio->lu , &ldA, ratio->ipiv, + ratio->Ai_y + i*np, &ldX, &info); + + /* Accumulate phase contributions */ + ratio->psum[i] = 0.0; + for (p = 0; p < np; p++) { + ratio->psum[i] += ratio->Ai_y[i*np + p]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_psys_contrib(grid_t *G, + well_t *W, + struct completion_data *wdata, + flowbc_t *bc, + struct compr_quantities *cq, + double dt, + const double *trans, + const double *gravcap_f, + const double *cpress0, + const double *porevol, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, nc, i, f; + size_t nconn; + double s; + + nc = G->number_of_cells; + nconn = G->cell_facepos[nc]; + + /* Compressible one-sided transmissibilities */ + set_dynamic_trans(G, trans, cq, h->pimpl->ratio); + compute_densrat_update(G, cq, h->pimpl->ratio, + h->pimpl->masstrans_f); + memcpy(h->pimpl->ctrans, + h->pimpl->ratio->psum, + nconn * sizeof *h->pimpl->ctrans); + + /* Compressible gravity contributions */ + set_dynamic_grav(G, bc, trans, gravcap_f, cq, h->pimpl->ratio); + compute_densrat_update(G, cq, h->pimpl->ratio, + h->pimpl->gravtrans_f); + + for (c = 0, i = 0; c < nc; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + s = 1.0 - 2.0*(G->face_cells[2*f + 0] != c); + + h->b[c] -= s * h->pimpl->ratio->psum[i]; + } + + h->b[c] += cq->voldiscr[c]; + } + + /* Compressible accumulation term (lhs and rhs) */ + compr_accum_term(nc, dt, porevol, cq->totcompr, h->pimpl->accum); + compr_src_add_press_accum(nc, cpress0, h->pimpl->accum, h->b); + + /* Compressible well contributions */ + if ((W != NULL) && (wdata != NULL)) { + nconn = W->well_connpos[ W->number_of_wells ]; + + set_dynamic_trans_well(W, cq->nphases, wdata, h->pimpl->ratio); + compute_densrat_update_well(W, wdata, cq, h->pimpl->ratio, + h->pimpl->masstrans_p); + memcpy(h->pimpl->wtrans, + h->pimpl->ratio->psum, nconn * sizeof *h->pimpl->wtrans); + + set_dynamic_grav_well(W, cq->nphases, wdata, h->pimpl->ratio); + compute_densrat_update_well(W, wdata, cq, h->pimpl->ratio, + h->pimpl->gravtrans_p); + memcpy(h->pimpl->wgpot, + h->pimpl->ratio->psum, nconn * sizeof *h->pimpl->wgpot); } } @@ -538,37 +628,107 @@ compute_cell_contrib(grid_t *G , /* ---------------------------------------------------------------------- */ static int assemble_cell_contrib(grid_t *G, - int c, + flowbc_t *bc, + const double *src, struct cfs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { - int c1, c2, i, f, j1, j2; + int c1, c2, c, i, f, j1, j2; int is_neumann; + const double *ctrans = h->pimpl->ctrans; + is_neumann = 1; - j1 = csrmatrix_elm_index(c, c, h->A); + for (c = i = 0; c < G->number_of_cells; c++) { + j1 = csrmatrix_elm_index(c, c, h->A); - h->A->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; + c2 = (c1 == c) ? c2 : c1; - c2 = (c1 == c) ? c2 : c1; + if (c2 >= 0) { + j2 = csrmatrix_elm_index(c, c2, h->A); - if (c2 >= 0) { - j2 = csrmatrix_elm_index(c, c2, h->A); + h->A->sa[j1] += ctrans[i]; + h->A->sa[j2] -= ctrans[i]; + } else if (bc->type[f] == PRESSURE) { + is_neumann = 0; - h->A->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; + h->A->sa[j1] += ctrans[i]; + h->b [c ] += ctrans[i] * bc->bcval[f]; + } + } + + h->b[c] += src[c]; + + /* Compressible accumulation term */ + h->A->sa[j1] += h->pimpl->accum[c]; + } + + return is_neumann; +} + + +/* ---------------------------------------------------------------------- */ +static int +assemble_well_contrib(size_t nc, + well_t *W, + well_control_t *wctrl, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, i, w; + int is_neumann, is_bhp; + + size_t jc, jw; + + double wtrans, dp; + + is_neumann = 1; + + for (w = i = 0; w < W->number_of_wells; w++) { + is_bhp = wctrl->ctrl[w] == BHP; + + for (; i < W->well_connpos[w + 1]; i++) { + c = W->well_cells[i]; + + wtrans = h->pimpl->wtrans[i]; /* WI * sum((Ac\Ap)*[pmob] */ + dp = h->pimpl->wgpot [i]; /* WI * sum((Ac\Ap)*[pmob]'*dP */ + + if (is_bhp) { + h->b[0 + c] += dp + wtrans * wctrl->target[w]; + h->b[nc + w] += wtrans * wctrl->target[w]; + } else { + jc = csrmatrix_elm_index(c, nc + w, h->A); + h->A->sa[jc] -= wtrans; + h->b [ c] += dp; + + jc = csrmatrix_elm_index(nc + w, c, h->A); + h->A->sa[jc] -= wtrans; + h->b[nc + w] -= dp; + } + + jc = csrmatrix_elm_index(0 + c, 0 + c, h->A); + jw = csrmatrix_elm_index(nc + w, nc + w, h->A); + + h->A->sa[jc] += wtrans; + h->A->sa[jw] += wtrans; + } + + is_neumann = is_neumann && (! is_bhp); + + if (! is_bhp) { + /* Enforce total (reservoir volume) rate constraint. */ + h->b[nc + w] += wctrl->target[w]; } } - h->b[ c ] = h->pimpl->ratio->residual; - - return 0; + return is_neumann; } @@ -678,37 +838,58 @@ compute_flux(grid_t *G, } -static size_t -maxconn(grid_t *G) +/* ---------------------------------------------------------------------- */ +static void +compute_wflux(well_t *W, + size_t np, + struct completion_data *wdata, + const double *cpress, + const double *wpress, + double *wflux) +/* ---------------------------------------------------------------------- */ { - int c; - size_t m, n; + int c, i, w; + size_t p; + double dp; - for (c = 0, m = 0; c < G->number_of_cells; c++) { - n = G->cell_facepos[c + 1] - G->cell_facepos[c]; + double *pmob, *gpot; - if (n > m) { m = n; } + pmob = wdata->phasemob; + gpot = wdata->gpot; + + for (w = i = 0; w < W->number_of_wells; w++) { + for (; i < W->well_connpos[w + 1]; i++) { + c = W->well_cells[i]; + + dp = wpress[w] - cpress[c]; + + wflux[i] = 0.0; + for (p = 0; p < np; p++) { + wflux[i] += pmob[i*np + p] * (dp + gpot[i*np + p]); + } + + wflux[i] *= wdata->WI[i]; + } } - - return m; } /* ---------------------------------------------------------------------- */ static int -is_incompr(int nc, const double *accum) +is_incompr(int nc, struct compr_quantities *cq) /* ---------------------------------------------------------------------- */ { int c, incompr; - for (c = 0, incompr = 1; (c < nc) && incompr; c++) { - incompr = ! (fabs(accum[c]) > 0.0); + for (c = 0, incompr = 1; (c < nc) && incompr; ++c) { + incompr = cq->totcompr[c] == 0.0; } return incompr; } + /* ====================================================================== * Public interface below separator. * ====================================================================== */ @@ -724,7 +905,7 @@ cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) new = malloc(1 * sizeof *new); if (new != NULL) { - new->pimpl = impl_allocate(G, W, maxconn(G), nphases); + new->pimpl = impl_allocate(G, W, nphases); new->A = construct_matrix(G, W); if ((new->pimpl == NULL) || (new->A == NULL)) { @@ -744,18 +925,23 @@ cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) } /* Allocate linear system components */ - new->b = new->pimpl->ddata + 0; - new->x = new->b + new->A->m; + new->b = new->pimpl->ddata + 0; + new->x = new->b + new->A->m; - new->pimpl->compflux_f = new->x + new->A->m; - new->pimpl->compflux_deriv_f = - new->pimpl->compflux_f + (nphases * nf); + /* Allocate reservoir components */ + new->pimpl->ctrans = new->x + new->A->m; + new->pimpl->accum = new->pimpl->ctrans + ngconn; + new->pimpl->masstrans_f = new->pimpl->accum + nc; + new->pimpl->gravtrans_f = new->pimpl->masstrans_f + (nphases * nf); - new->pimpl->flux_work = - new->pimpl->compflux_deriv_f + (nphases * 2 * nf); + /* Allocate well components */ + new->pimpl->wtrans = new->pimpl->gravtrans_f + (nphases * nf); + new->pimpl->wgpot = new->pimpl->wtrans + nwconn; + new->pimpl->masstrans_p = new->pimpl->wgpot + nwconn; + new->pimpl->gravtrans_p = new->pimpl->masstrans_p + (nphases * nwconn); - new->pimpl->scratch_f = - new->pimpl->flux_work + (nphases * (1 + 2)); + /* Allocate scratch array for face pressure calculations */ + new->pimpl->scratch_f = new->pimpl->gravtrans_p + (nphases * nwconn); } return new; @@ -766,39 +952,27 @@ cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) void cfs_tpfa_assemble(grid_t *G, double dt, + well_t *W, flowbc_t *bc, const double *src, - const double *zc, struct compr_quantities *cq, const double *trans, const double *gravcap_f, - const double *cpress, + well_control_t *wctrl, + struct completion_data *wdata, + const double *cpress0, const double *porevol, struct cfs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { - int res_is_neumann, c, np2; + int res_is_neumann, well_is_neumann; csrmatrix_zero( h->A); vector_zero (h->A->m, h->b); - h->pimpl->is_incomp = 1; + compute_psys_contrib(G, W, wdata, bc, cq, dt, + trans, gravcap_f, cpress0, porevol, h); - compute_compflux_and_deriv(G, cq->nphases, cpress, trans, gravcap_f, - cq->phasemobf, cq->Af, h->pimpl); - - np2 = cq->nphases * cq->nphases; - for (c = 0; c < G->number_of_cells; - c++, zc += cq->nphases) { - - compute_cell_contrib(G, c, cq->nphases, porevol[c], dt, zc, - cq->Ac + (c * np2), cq->dAc + (c * np2), - h->pimpl); - - assemble_cell_contrib(G, c, h); - } - -#if 0 res_is_neumann = assemble_cell_contrib(G, bc, src, h); if ((W != NULL) && (wctrl != NULL)) { @@ -808,18 +982,14 @@ cfs_tpfa_assemble(grid_t *G, } else { well_is_neumann = 1; } -#endif - res_is_neumann = 1; - - if (res_is_neumann && h->pimpl->is_incomp) { - /*&& well_is_neumann && h->pimpl->is_incomp) {*/ + if (res_is_neumann && well_is_neumann && + is_incompr(G->number_of_cells, cq)) { h->A->sa[0] *= 2; } } -#if 0 /* ---------------------------------------------------------------------- */ void cfs_tpfa_press_increment(grid_t *G, @@ -1129,7 +1299,6 @@ cfs_tpfa_expl_mass_transport(grid_t *G, } } } -#endif /* ---------------------------------------------------------------------- */ diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index e96c1e938..7cbec29a6 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -47,32 +47,38 @@ cfs_tpfa_construct(grid_t *G, well_t *W, int nphases); void cfs_tpfa_assemble(grid_t *G, double dt, + well_t *W, flowbc_t *bc, const double *src, - const double *zc, struct compr_quantities *cq, const double *trans, const double *gravcap_f, - const double *cpress, + well_control_t *wctrl, + struct completion_data *wdata, + const double *cpress0, const double *porevol, struct cfs_tpfa_data *h); void cfs_tpfa_press_increment(grid_t *G, + well_t *W, struct cfs_tpfa_data *h, - double *cpress_inc); + double *cpress_inc, + double *wpress_inc); -#if 0 void cfs_tpfa_flux(grid_t *G, flowbc_t *bc, + well_t *W, int np, const double *trans, const double *pmobf, const double *gravcap_f, const double *cpress, - double *fflux); -#endif + const double *wpress, + struct completion_data *wdata, + double *fflux, + double *wflux); void cfs_tpfa_fpress(grid_t *G,