From 40debf6a5fb56cc690da07830da8df0ddbc6d68d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 7 Dec 2010 17:12:09 +0100 Subject: [PATCH] Add tentative Peaceman well support. Callers may pass NULLs in absence of wells in any given model. This implementation assembles an equation for each well, irrespective of well control type (BHP or RATE), and assumes that productivity indices and perforation pressure drops account for multiphase effects. --- src/cfs_tpfa.c | 147 +++++++++++++++++++++++++++++++++++++------------ src/cfs_tpfa.h | 7 ++- 2 files changed, 118 insertions(+), 36 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 4d8cb7d38..194aa13a0 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -458,6 +458,106 @@ compute_psys_contrib(grid_t *G, } +/* ---------------------------------------------------------------------- */ +static int +assemble_cell_contrib(grid_t *G, + flowbc_t *bc, + const double *src, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2, c, i, f, j1, j2; + int is_neumann; + + const double *ctrans = h->pimpl->dd->ctrans; + + is_neumann = 0; + + for (c = i = 0; c < G->number_of_cells; c++) { + j1 = csrmatrix_elm_index(c, c, h->A); + + 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]; + + c2 = (c1 == c) ? c2 : c1; + + 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[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->dd->P[c]; + } + + return is_neumann; +} + + +/* ---------------------------------------------------------------------- */ +static int +assemble_well_contrib(size_t nc, + well_t *W, + well_control_t *wctrl, + const double *WI, + const double *wdp, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, i, w; + int is_neumann, is_bhp, is_inj; + + size_t jc, jw; + + is_neumann = 1; + + for (w = i = 0; w < W->number_of_wells; w++) { + is_bhp = wctrl->ctrl[w] == BHP; + is_inj = wctrl->type[w] == INJECTOR; + + for (; i < W->well_connpos[w + 1]; i++) { + c = W->well_cells[i]; + + if (is_bhp) { + h->b[0 + c] += WI[i] * (wctrl->target[w] + wdp[i]); + h->b[nc + w] += WI[i] * wctrl->target[w]; + } else { + jc = csrmatrix_elm_index(c, nc + w, h->A); + h->A->sa[jc] -= WI[i]; + h->b [ c] += WI[i] * wdp[i]; + + jc = csrmatrix_elm_index(nc + w, c, h->A); + h->A->sa[jc] -= WI[i]; + h->b[nc + w] -= WI[i] * wdp[i]; + } + + jc = csrmatrix_elm_index(0 + c, 0 + c, h->A); + jw = csrmatrix_elm_index(nc + w, nc + w, h->A); + + h->A->sa[jc] += WI[i]; + h->A->sa[jw] += WI[i]; + } + + is_neumann = is_neumann && (! is_bhp); + } + + return is_neumann; +} + + /* ---------------------------------------------------------------------- */ static void compute_fpress(grid_t *G, @@ -560,7 +660,7 @@ compute_flux(grid_t *G, /* ---------------------------------------------------------------------- */ struct cfs_tpfa_data * -cfs_tpfa_construct(grid_t *G, int nphases) +cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) /* ---------------------------------------------------------------------- */ { size_t nf; @@ -569,8 +669,8 @@ cfs_tpfa_construct(grid_t *G, int nphases) new = malloc(1 * sizeof *new); if (new != NULL) { - new->pimpl = impl_allocate(G, NULL, nphases); - new->A = construct_matrix(G, NULL); + new->pimpl = impl_allocate(G, W, nphases); + new->A = construct_matrix(G, W); if ((new->pimpl == NULL) || (new->A == NULL)) { cfs_tpfa_destroy(new); @@ -598,56 +698,33 @@ cfs_tpfa_construct(grid_t *G, int nphases) void cfs_tpfa_assemble(grid_t *G, double dt, + well_t *W, flowbc_t *bc, const double *src, struct compr_quantities *cq, const double *trans, const double *gravcap_f, + well_control_t *wctrl, + const double *WI, + const double *wdp, const double *cpress0, const double *porevol, struct cfs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { - int c1, c2, c, i, f, j1, j2; int is_neumann; - double *ctrans = h->pimpl->dd->ctrans; - csrmatrix_zero( h->A); vector_zero (h->A->m, h->b); compute_psys_contrib(G, cq, dt, trans, gravcap_f, cpress0, porevol, h); - is_neumann = 1; + is_neumann = assemble_cell_contrib(G, bc, src, h); - for (c = i = 0; c < G->number_of_cells; c++) { - j1 = csrmatrix_elm_index(c, c, h->A); - - 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]; - - c2 = (c1 == c) ? c2 : c1; - - 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[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->dd->P[c]; + if ((W != NULL) && (wctrl != NULL) && + (WI != NULL) && (wdp != NULL)) { + is_neumann = is_neumann && + assemble_well_contrib(G->number_of_cells, W, wctrl, WI, wdp, h); } if (is_neumann) { diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index 9b9cff2bc..976333325 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -22,6 +22,7 @@ #include "grid.h" #include "flow_bc.h" +#include "well.h" #ifdef __cplusplus extern "C" { @@ -41,16 +42,20 @@ struct cfs_tpfa_data { struct cfs_tpfa_data * -cfs_tpfa_construct(grid_t *G, int nphases); +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, struct compr_quantities *cq, const double *trans, const double *gravcap_f, + well_control_t *wctrl, + const double *WI, + const double *wdp, const double *cpress0, const double *porevol, struct cfs_tpfa_data *h);