From 9b0b2581980667ff77f97d674e6a24e3e74f835c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 17 Oct 2011 11:05:04 +0200 Subject: [PATCH 01/36] Factor pressure (increment) assignment out of _press_flux(). The linear solution h->x is the pressure increment, not the actual pressure value, so we cannot compute fluxes based on h->x alone. --- src/cfs_tpfa.c | 50 ++++++++++++++++++++++++++++++++------------------ src/cfs_tpfa.h | 32 +++++++++++++++++++------------- 2 files changed, 51 insertions(+), 31 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 90d58eeb..0fea19f9 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -992,34 +992,48 @@ cfs_tpfa_assemble(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_press_flux(grid_t *G, - flowbc_t *bc, - well_t *W, - int np, - const double *trans, - const double *pmobf, - const double *gravcap_f, - struct completion_data *wdata, - struct cfs_tpfa_data *h, - double *cpress, - double *fflux, - double *wpress, - double *wflux) +cfs_tpfa_press_increment(grid_t *G, + well_t *W, + struct cfs_tpfa_data *h, + double *cpress_inc, + double *wpress_inc) /* ---------------------------------------------------------------------- */ { /* Assign cell pressure directly from solution vector */ - memcpy(cpress, h->x, G->number_of_cells * sizeof *cpress); + memcpy(cpress_inc, h->x, G->number_of_cells * sizeof *cpress_inc); + if (W != NULL) { + assert (wpress_inc != NULL); + + /* Assign well BHP directly from solution vector */ + memcpy(wpress_inc, h->x + G->number_of_cells, + W->number_of_wells * sizeof *wpress_inc); + } +} + + +/* ---------------------------------------------------------------------- */ +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, + const double *wpress, + struct completion_data *wdata, + double *fflux, + double *wflux) +/* ---------------------------------------------------------------------- */ +{ compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux); if ((W != NULL) && (wdata != NULL)) { assert (wpress != NULL); assert (wflux != NULL); - /* Assign well BHP directly from solution vector */ - memcpy(wpress, h->x + G->number_of_cells, - W->number_of_wells * sizeof *wpress); - compute_wflux(W, np, wdata, cpress, wpress, wflux); } } diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index acd809ea..7cbec29a 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -60,19 +60,25 @@ cfs_tpfa_assemble(grid_t *G, struct cfs_tpfa_data *h); void -cfs_tpfa_press_flux(grid_t *G, - flowbc_t *bc, - well_t *W, - int np, - const double *trans, - const double *pmobf, - const double *gravcap_f, - struct completion_data *wdata, - struct cfs_tpfa_data *h, - double *cpress, - double *fflux, - double *wpress, - double *wflux); +cfs_tpfa_press_increment(grid_t *G, + well_t *W, + struct cfs_tpfa_data *h, + double *cpress_inc, + double *wpress_inc); + +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, + const double *wpress, + struct completion_data *wdata, + double *fflux, + double *wflux); void cfs_tpfa_fpress(grid_t *G, From ceaea9f98700f6324254c0b49b6b118c27b1331d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:38:28 +0200 Subject: [PATCH 02/36] Rework compressibility representation. Switch to storing a complete fluid-matrix derivative in the compr_quantities rather than the total compressibility. Maintain the "volume discrepancy" field. Also, add traditional memory management functions. --- src/compr_quant.c | 52 +++++++++++++++++++++++++++++++++++++++++++++++ src/compr_quant.h | 12 ++++++++--- 2 files changed, 61 insertions(+), 3 deletions(-) diff --git a/src/compr_quant.c b/src/compr_quant.c index 296168c7..93d12a75 100644 --- a/src/compr_quant.c +++ b/src/compr_quant.c @@ -17,9 +17,61 @@ along with OPM. If not, see . */ +#include + +#include "sparse_sys.h" #include "compr_quant.h" +void +compr_quantities_deallocate(struct compr_quantities *cq) +{ + if (cq != NULL) { + free(cq->Ac); + } + + free(cq); +} + + +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np) +{ + size_t alloc_sz, np2; + struct compr_quantities *cq; + + cq = malloc(1 * sizeof *cq); + + if (cq != NULL) { + np2 = np * np; + + alloc_sz = np2 * nc; /* Ac */ + alloc_sz += np2 * nc; /* dAc */ + alloc_sz += np2 * nf; /* Af */ + alloc_sz += np * nf; /* phasemobf */ + alloc_sz += 1 * nc; /* voldiscr */ + + cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); + + if (cq->Ac == NULL) { + compr_quantities_deallocate(cq); + cq = NULL; + } else { + cq->dAc = cq->Ac + (np2 * nc); + cq->Af = cq->dAc + (np2 * nc); + cq->phasemobf = cq->Af + (np2 * nf); + cq->voldiscr = cq->phasemobf + (np * nf); + + cq->nphases = np; + + vector_zero(alloc_sz, cq->Ac); + } + } + + return cq; +} + + /* ---------------------------------------------------------------------- */ /* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ /* ---------------------------------------------------------------------- */ diff --git a/src/compr_quant.h b/src/compr_quant.h index cacfb551..b3531176 100644 --- a/src/compr_quant.h +++ b/src/compr_quant.h @@ -29,15 +29,21 @@ extern "C" { #endif struct compr_quantities { - int nphases; /* Number of phases/components */ + int nphases; /* Number of phases/components */ - double *totcompr; /* Total compressibility per cell */ - double *voldiscr; /* Volume discrepancy per cell */ double *Ac; /* RB^{-1} per cell */ + double *dAc; /* d/dp (RB^{-1}) per cell */ double *Af; /* RB^{-1} per face */ double *phasemobf; /* Phase mobility per face */ + double *voldiscr; /* Volume discrepancy per cell */ }; +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np); + +void +compr_quantities_deallocate(struct compr_quantities *cq); + void compr_flux_term(grid_t *G, const double *fflux, From e44b7317b8473eb26e334df2c6f090798aa44135 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:40:19 +0200 Subject: [PATCH 03/36] Delete trailing whitespace. --- src/sparse_sys.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sparse_sys.c b/src/sparse_sys.c index d3e97df6..793941bf 100644 --- a/src/sparse_sys.c +++ b/src/sparse_sys.c @@ -247,7 +247,7 @@ vector_write(size_t n, const double *v, const char *fn) if (fp != NULL) { vector_write_stream(n, v, fp); } - + fclose(fp); } From 99fa7b7f0ed5d0e67ea3bb21c91414fef69af982 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:41:13 +0200 Subject: [PATCH 04/36] Delete trailing whitespace. --- src/blas_lapack.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blas_lapack.h b/src/blas_lapack.h index 4c0230b6..7051783a 100644 --- a/src/blas_lapack.h +++ b/src/blas_lapack.h @@ -70,7 +70,7 @@ void dgetrs_(const char *trans, const MAT_SIZE_T *n, const double *A , const MAT_SIZE_T *lda, const MAT_SIZE_T *ipiv , double *B, const MAT_SIZE_T *ldb , MAT_SIZE_T *info); - + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From cd02647d53fc8889e54dfb111b7b5f7403545f2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:41:44 +0200 Subject: [PATCH 05/36] Delete trailing whitespace. --- src/test_cfs_tpfa.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/test_cfs_tpfa.c b/src/test_cfs_tpfa.c index 6becf0ae..ab59c11f 100644 --- a/src/test_cfs_tpfa.c +++ b/src/test_cfs_tpfa.c @@ -674,7 +674,7 @@ main(void) for (i = 0; i < G->number_of_cells; i++) { fprintf(stderr, "press(%02d) = %g;\n", i + 1, cpress[i]); } - + cfs_tpfa_destroy(h); deallocate_cq(cq); deallocate_flowbc(bc); From 4f657d9f872c96e036f484462446bc9440ed8420 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:46:05 +0200 Subject: [PATCH 06/36] Make first attempt at implementing the residual/Jacobian formulation. This is a step in the direction of having a true IMPES pressure solver and to remove the 'experimental_jacobian' option in dune-porsol. --- src/cfs_tpfa.c | 813 ++++++++++++++++++++----------------------------- src/cfs_tpfa.h | 18 +- 2 files changed, 328 insertions(+), 503 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 0fea19f9..8ce3efc0 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -24,33 +24,25 @@ struct densrat_util { MAT_SIZE_T *ipiv; + + double residual; double *lu; - - double *x; - double *Ai_y; - double *psum; - - /* Storage */ - double *ddata; + double *t1; + double *t2; + double *mat_row; + double *coeff; + double *linsolve_buffer; }; struct cfs_tpfa_impl { - /* Reservoir flow */ - double *ctrans; - double *accum; + int is_incomp; /* One entry per component per face */ - double *masstrans_f; /* RB^{-1} [ phase-mobility ] */ - double *gravtrans_f; /* RB^{-1} [ grav + capillary ] */ + double *compflux_f; /* A_{ij} v_{ij} */ + double *compflux_deriv_f; /* A_{ij} \partial_{p} v_{ij} */ - /* 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 ] */ + double *flux_work; /* Scratch array for face pressure calculation */ double *scratch_f; @@ -68,7 +60,7 @@ deallocate_densrat(struct densrat_util *ratio) /* ---------------------------------------------------------------------- */ { if (ratio != NULL) { - free(ratio->ddata); + free(ratio->lu); free(ratio->ipiv); } @@ -78,46 +70,41 @@ deallocate_densrat(struct densrat_util *ratio) /* ---------------------------------------------------------------------- */ static struct densrat_util * -allocate_densrat(grid_t *g, well_t *w, int np) +allocate_densrat(size_t max_conn, int np) /* ---------------------------------------------------------------------- */ { - int ntotperf; - size_t nglobconn, ntotconn, ddata_sz; + size_t alloc_sz, n_buffer_col; + struct densrat_util *ratio; - struct densrat_util *new; + ratio = malloc(1 * sizeof *ratio); - new = malloc(1 * sizeof *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} */ - if (new != NULL) { - if (w != NULL) { - ntotperf = w->well_connpos[ w->number_of_wells ]; + 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; } else { - 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 ; + 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 ); } } - return new; + return ratio; } @@ -137,7 +124,7 @@ impl_deallocate(struct cfs_tpfa_impl *pimpl) /* ---------------------------------------------------------------------- */ static struct cfs_tpfa_impl * -impl_allocate(grid_t *G, well_t *W, int np) +impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) /* ---------------------------------------------------------------------- */ { size_t nnu, ngconn, nwperf; @@ -158,24 +145,18 @@ impl_allocate(grid_t *G, well_t *W, int np) ddata_sz = 2 * nnu; /* b, x */ /* Reservoir */ - 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 * G->number_of_faces ; /* compflux_f */ + ddata_sz += np * (2 * G->number_of_faces); /* compflux_deriv_f */ - /* 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 += np * (1 + 2) ; /* flux_work */ - 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(G, W, np); + new->ratio = allocate_densrat(max_conn, np); if (new->ddata == NULL || new->ratio == NULL) { impl_deallocate(new); @@ -279,348 +260,277 @@ construct_matrix(grid_t *G, well_t *W) } -/* ---------------------------------------------------------------------- */ static void -solve_cellsys_core(grid_t *G , - size_t sz , - const double *Ac , - const double *bf , - double *xcf , - double *luAc, - MAT_SIZE_T *ipiv) -/* ---------------------------------------------------------------------- */ +factorise_fluid_matrix(int np, const double *A, struct densrat_util *ratio) { - int c, i, f; - size_t j, p2; - double *v; + int np2; + MAT_SIZE_T m, n, ld, info; - MAT_SIZE_T nrows, ncols, ldA, ldX, nrhs, info; + m = n = ld = np; + np2 = np * np; - nrows = ncols = ldA = ldX = sz; - info = 0; + memcpy (ratio->lu, A, np2 * sizeof *ratio->lu); + dgetrf_(&m, &n, ratio->lu, &ld, ratio->ipiv, &info); - 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; - } + assert (info == 0); } -/* ---------------------------------------------------------------------- */ static void -small_matvec(size_t n, - int sz, - const double *A, - const double *X, - double *Y) -/* ---------------------------------------------------------------------- */ +solve_linear_systems(int np , + MAT_SIZE_T nrhs , + struct densrat_util *ratio, + double *b ) { - size_t i, p1, p2; + MAT_SIZE_T n, ldA, ldB, info; - MAT_SIZE_T nrows, ncols, ld, incx, incy; - double a1, a2; + n = ldA = ldB = np; - nrows = ncols = ld = sz; - incx = incy = 1; + dgetrs_("No Transpose", &n, + &nrhs, ratio->lu, &ldA, ratio->ipiv, + b , &ldB, &info); + 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; - for (i = p1 = p2 = 0; i < n; i++) { - dgemv_("No Transpose", &nrows, &ncols, - &a1, A + p2, &ld, X + p1, &incx, - &a2, Y + p1, &incy); + dgemm_("No Transpose", "No Transpose", &m, &n, &k, + &a1, A, &ldA, B, &ldB, &a2, C, &ldC); +} - 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) +{ + 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; } } -/* ---------------------------------------------------------------------- */ static void -solve_cellsys(grid_t *G , - size_t sz, - const double *Ac, - const double *bf, - struct densrat_util *ratio) -/* ---------------------------------------------------------------------- */ +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 ) { - solve_cellsys_core(G, sz, Ac, bf, ratio->Ai_y, - ratio->lu, ratio->ipiv); -} + int c1, c2, f, np2; + double dp; + double *cflux, *dcflux; + np2 = np * np; -/* ---------------------------------------------------------------------- */ -static void -set_dynamic_trans(grid_t *G , - const double *trans, - struct compr_quantities *cq , - struct densrat_util *ratio) -/* ---------------------------------------------------------------------- */ -{ - int f, p, i; + cflux = pimpl->compflux_f; + dcflux = pimpl->compflux_deriv_f; - 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]; - } - } -} + for (f = 0; f < G->number_of_faces; + f++, pmobf += np, gcapf += np, Af += np2, + cflux += np, dcflux += 2 * np) { - -/* ---------------------------------------------------------------------- */ -static void -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 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; - } + 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 void -set_dynamic_trans_well(well_t *W, - size_t np, - struct completion_data *wdata, - struct densrat_util *ratio) -/* ---------------------------------------------------------------------- */ +static int +count_internal_conn(grid_t *G, int c) { - size_t p, i, k, nconn; + int c1, c2, f, i, nconn; - nconn = W->well_connpos[ W->number_of_wells ]; + 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]; - for (i = k = 0; i < nconn; i++) { - for (p = 0; p < np; p++, k++) { - ratio->x[k] = wdata->WI[i] * wdata->phasemob[k]; + 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; } } + + assert (conn == nconn + 1); + assert (cflx == pimpl->ratio->linsolve_buffer + (nconn + 1)*np); + + return nconn; } -/* ---------------------------------------------------------------------- */ static void -set_dynamic_grav_well(well_t *W, - size_t np, - struct completion_data *wdata, - struct densrat_util *ratio) -/* ---------------------------------------------------------------------- */ +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) { - size_t p, i, k, nconn; + int c1, c2, f, i, nconn, np2, p; + MAT_SIZE_T nrhs; + double s, dF1, dF2, *dv, *dv1, *dv2; - 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; - nrows = ncols = ldA = ldX = np; - incx = incy = nrhs = 1; + nconn = init_cell_contrib(G, c, np, pvol, dt, z, pimpl); + nrhs = 1 + (1 + 2)*nconn; /* [z, Af*v, Af*dv] */ - a1 = 1.0; - a2 = 0.0; + factorise_fluid_matrix(np, Ac, pimpl->ratio); + solve_linear_systems (np, nrhs, pimpl->ratio, + pimpl->ratio->linsolve_buffer); - for (i = 0; i < nconn; i++) { - c = W->well_cells[i]; + /* 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); - /* 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); - - /* 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]; + /* Compute residual in cell 'c' */ + pimpl->ratio->residual = pvol; + for (p = 0; p < np; p++) { + pimpl->ratio->residual += pimpl->ratio->t1[ p ]; } - /* 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); + /* Jacobian row */ - /* Compressible well contributions */ - if ((W != NULL) && (wdata != NULL)) { - nconn = W->well_connpos[ W->number_of_wells ]; + /* 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); - 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); + dF1 = dF2 = 0.0; + for (p = 0; p < np; p++) { + dF1 += pimpl->ratio->t1[ p ]; + dF2 += pimpl->ratio->t2[ p ]; + } - 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); + 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; + } } } @@ -628,107 +538,37 @@ compute_psys_contrib(grid_t *G, /* ---------------------------------------------------------------------- */ static int assemble_cell_contrib(grid_t *G, - flowbc_t *bc, - const double *src, + int c, struct cfs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { - int c1, c2, c, i, f, j1, j2; + int c1, c2, i, f, j1, j2; int is_neumann; - const double *ctrans = h->pimpl->ctrans; - is_neumann = 1; - for (c = i = 0; c < G->number_of_cells; c++) { - j1 = csrmatrix_elm_index(c, c, h->A); + j1 = csrmatrix_elm_index(c, c, h->A); - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; + h->A->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; - c2 = (c1 == c) ? c2 : c1; + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; - if (c2 >= 0) { - j2 = csrmatrix_elm_index(c, c2, h->A); + c2 = (c1 == c) ? c2 : c1; - h->A->sa[j1] += ctrans[i]; - h->A->sa[j2] -= ctrans[i]; - } else if (bc->type[f] == PRESSURE) { - is_neumann = 0; + if (c2 >= 0) { + j2 = csrmatrix_elm_index(c, c2, h->A); - 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->A->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; } } - return is_neumann; + h->b[ c ] = h->pimpl->ratio->residual; + + return 0; } @@ -838,58 +678,37 @@ compute_flux(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) -/* ---------------------------------------------------------------------- */ +static size_t +maxconn(grid_t *G) { - int c, i, w; - size_t p; - double dp; + int c; + size_t m, n; - double *pmob, *gpot; + for (c = 0, m = 0; c < G->number_of_cells; c++) { + n = G->cell_facepos[c + 1] - G->cell_facepos[c]; - 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]; - } + if (n > m) { m = n; } } + + return m; } /* ---------------------------------------------------------------------- */ static int -is_incompr(int nc, struct compr_quantities *cq) +is_incompr(int nc, const double *accum) /* ---------------------------------------------------------------------- */ { int c, incompr; - for (c = 0, incompr = 1; (c < nc) && incompr; ++c) { - incompr = cq->totcompr[c] == 0.0; + for (c = 0, incompr = 1; (c < nc) && incompr; c++) { + incompr = ! (fabs(accum[c]) > 0.0); } return incompr; } - /* ====================================================================== * Public interface below separator. * ====================================================================== */ @@ -905,7 +724,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, nphases); + new->pimpl = impl_allocate(G, W, maxconn(G), nphases); new->A = construct_matrix(G, W); if ((new->pimpl == NULL) || (new->A == NULL)) { @@ -925,23 +744,18 @@ 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; - /* 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->compflux_f = new->x + new->A->m; + new->pimpl->compflux_deriv_f = + new->pimpl->compflux_f + (nphases * 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->flux_work = + new->pimpl->compflux_deriv_f + (nphases * 2 * nf); - /* Allocate scratch array for face pressure calculations */ - new->pimpl->scratch_f = new->pimpl->gravtrans_p + (nphases * nwconn); + new->pimpl->scratch_f = + new->pimpl->flux_work + (nphases * (1 + 2)); } return new; @@ -952,27 +766,39 @@ 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, - well_control_t *wctrl, - struct completion_data *wdata, - const double *cpress0, + const double *cpress, const double *porevol, struct cfs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { - int res_is_neumann, well_is_neumann; + int res_is_neumann, c, np2; csrmatrix_zero( h->A); vector_zero (h->A->m, h->b); - compute_psys_contrib(G, W, wdata, bc, cq, dt, - trans, gravcap_f, cpress0, porevol, h); + h->pimpl->is_incomp = 1; + 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)) { @@ -982,14 +808,18 @@ cfs_tpfa_assemble(grid_t *G, } else { well_is_neumann = 1; } +#endif - if (res_is_neumann && well_is_neumann && - is_incompr(G->number_of_cells, cq)) { + res_is_neumann = 1; + + if (res_is_neumann && h->pimpl->is_incomp) { + /*&& well_is_neumann && h->pimpl->is_incomp) {*/ h->A->sa[0] *= 2; } } +#if 0 /* ---------------------------------------------------------------------- */ void cfs_tpfa_press_increment(grid_t *G, @@ -1299,6 +1129,7 @@ cfs_tpfa_expl_mass_transport(grid_t *G, } } } +#endif /* ---------------------------------------------------------------------- */ diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index 7cbec29a..e96c1e93 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -47,38 +47,32 @@ 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, - well_control_t *wctrl, - struct completion_data *wdata, - const double *cpress0, + const double *cpress, 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 *wpress_inc); + double *cpress_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, - const double *wpress, - struct completion_data *wdata, - double *fflux, - double *wflux); + double *fflux); +#endif void cfs_tpfa_fpress(grid_t *G, From 7bea27208e564c5e2126d0be933092fd74dec6e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:12:54 +0200 Subject: [PATCH 07/36] Grab copy of 'cfs_tpfa' module in preparation of restoring backwards compat. --- src/cfs_tpfa_residual.c | 1146 +++++++++++++++++++++++++++++++++++++++ src/cfs_tpfa_residual.h | 128 +++++ 2 files changed, 1274 insertions(+) create mode 100644 src/cfs_tpfa_residual.c create mode 100644 src/cfs_tpfa_residual.h diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c new file mode 100644 index 00000000..8ce3efc0 --- /dev/null +++ b/src/cfs_tpfa_residual.c @@ -0,0 +1,1146 @@ +#include +#include +#include +#include +#include + +#include "blas_lapack.h" +#include "flow_bc.h" +#include "well.h" + +#include "compr_quant.h" +#include "trans_tpfa.h" +#include "cfs_tpfa.h" +#include "sparse_sys.h" + + +#if defined(MAX) +#undef MAX +#endif + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + + + +struct densrat_util { + MAT_SIZE_T *ipiv; + + double residual; + double *lu; + double *t1; + double *t2; + double *mat_row; + double *coeff; + double *linsolve_buffer; +}; + + +struct cfs_tpfa_impl { + int is_incomp; + + /* One entry per component per face */ + double *compflux_f; /* A_{ij} v_{ij} */ + double *compflux_deriv_f; /* A_{ij} \partial_{p} v_{ij} */ + + double *flux_work; + + /* Scratch array for face pressure calculation */ + double *scratch_f; + + struct densrat_util *ratio; + + /* Linear storage */ + double *ddata; +}; + + +/* ---------------------------------------------------------------------- */ +static void +deallocate_densrat(struct densrat_util *ratio) +/* ---------------------------------------------------------------------- */ +{ + if (ratio != NULL) { + free(ratio->lu); + free(ratio->ipiv); + } + + free(ratio); +} + + +/* ---------------------------------------------------------------------- */ +static struct densrat_util * +allocate_densrat(size_t max_conn, int np) +/* ---------------------------------------------------------------------- */ +{ + size_t alloc_sz, n_buffer_col; + struct densrat_util *ratio; + + ratio = malloc(1 * sizeof *ratio); + + 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} */ + + 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; + } 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 ); + } + } + + return ratio; +} + + +/* ---------------------------------------------------------------------- */ +static void +impl_deallocate(struct cfs_tpfa_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + if (pimpl != NULL) { + free (pimpl->ddata); + deallocate_densrat(pimpl->ratio); + } + + free(pimpl); +} + + +/* ---------------------------------------------------------------------- */ +static struct cfs_tpfa_impl * +impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) +/* ---------------------------------------------------------------------- */ +{ + size_t nnu, ngconn, nwperf; + struct cfs_tpfa_impl *new; + + size_t ddata_sz; + + nnu = G->number_of_cells; + ngconn = G->cell_facepos[ G->number_of_cells ]; + nwperf = 0; + + if (W != NULL) { + nnu += W->number_of_wells; + nwperf = W->well_connpos[ W->number_of_wells ]; + } + + /* Linear system */ + 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 += np * (1 + 2) ; /* flux_work */ + + 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); + + if (new->ddata == NULL || new->ratio == NULL) { + impl_deallocate(new); + new = NULL; + } + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +static struct CSRMatrix * +construct_matrix(grid_t *G, well_t *W) +/* ---------------------------------------------------------------------- */ +{ + int f, c1, c2, w, i, nc, nnu; + size_t nnz; + + struct CSRMatrix *A; + + nc = nnu = G->number_of_cells; + if (W != NULL) { + nnu += W->number_of_wells; + } + + A = csrmatrix_new_count_nnz(nnu); + + if (A != NULL) { + /* Self connections */ + for (i = 0; i < nnu; i++) { + A->ia[ i + 1 ] = 1; + } + + /* Other connections */ + for (f = 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)) { + A->ia[ c1 + 1 ] += 1; + A->ia[ c2 + 1 ] += 1; + } + } + + if (W != NULL) { + /* Well <-> cell connections */ + for (w = i = 0; w < W->number_of_wells; w++) { + for (; i < W->well_connpos[w + 1]; i++) { + c1 = W->well_cells[i]; + + A->ia[ 0 + c1 + 1 ] += 1; /* c -> w */ + A->ia[ nc + w + 1 ] += 1; /* w -> c */ + } + } + } + + nnz = csrmatrix_new_elms_pushback(A); + if (nnz == 0) { + csrmatrix_delete(A); + A = NULL; + } + } + + if (A != NULL) { + /* Fill self connections */ + for (i = 0; i < nnu; i++) { + A->ja[ A->ia[ i + 1 ] ++ ] = i; + } + + /* Fill other connections */ + for (f = 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)) { + A->ja[ A->ia[ c1 + 1 ] ++ ] = c2; + A->ja[ A->ia[ c2 + 1 ] ++ ] = c1; + } + } + + if (W != NULL) { + /* Fill well <-> cell connections */ + for (w = i = 0; w < W->number_of_wells; w++) { + for (; i < W->well_connpos[w + 1]; i++) { + c1 = W->well_cells[i]; + + A->ja[ A->ia[ 0 + c1 + 1 ] ++ ] = nc + w; + A->ja[ A->ia[ nc + w + 1 ] ++ ] = c1 ; + } + } + } + + assert ((size_t) A->ia[ nnu ] == nnz); + + /* Enforce sorted connection structure per row */ + csrmatrix_sortrows(A); + } + + return A; +} + + +static void +factorise_fluid_matrix(int np, const double *A, struct densrat_util *ratio) +{ + int np2; + MAT_SIZE_T m, n, ld, info; + + m = n = ld = np; + np2 = np * np; + + memcpy (ratio->lu, A, np2 * sizeof *ratio->lu); + dgetrf_(&m, &n, ratio->lu, &ld, ratio->ipiv, &info); + + assert (info == 0); +} + + +static void +solve_linear_systems(int np , + MAT_SIZE_T nrhs , + struct densrat_util *ratio, + double *b ) +{ + MAT_SIZE_T n, ldA, ldB, info; + + n = ldA = ldB = np; + + dgetrs_("No Transpose", &n, + &nrhs, ratio->lu, &ldA, ratio->ipiv, + b , &ldB, &info); + + 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); +} + + +static void +compute_darcyflux_and_deriv(int np, + double trans, + double dp, + const double *pmobf, + const double *gcapf, + double *dflux, + double *dflux_deriv) +{ + 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; + } +} + + +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 ) +{ + int c1, c2, f, np2; + double dp; + double *cflux, *dcflux; + + 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; + } + } + + 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) +{ + int c1, c2, f, i, nconn, np2, p; + MAT_SIZE_T nrhs; + double s, dF1, dF2, *dv, *dv1, *dv2; + + np2 = np * np; + + nconn = init_cell_contrib(G, c, np, pvol, dt, z, pimpl); + nrhs = 1 + (1 + 2)*nconn; /* [z, Af*v, Af*dv] */ + + factorise_fluid_matrix(np, Ac, pimpl->ratio); + solve_linear_systems (np, nrhs, pimpl->ratio, + pimpl->ratio->linsolve_buffer); + + /* 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); + + /* Compute residual in cell 'c' */ + pimpl->ratio->residual = pvol; + for (p = 0; p < np; p++) { + pimpl->ratio->residual += pimpl->ratio->t1[ p ]; + } + + /* 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; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static int +assemble_cell_contrib(grid_t *G, + int c, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2, i, f, j1, j2; + int is_neumann; + + is_neumann = 1; + + 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]; + + 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[j2] -= h->pimpl->ratio->mat_row[ i ]; + } + } + + h->b[ c ] = h->pimpl->ratio->residual; + + return 0; +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + const double *cpress, + const double *fflux, + double *fpress, + double *scratch_f) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f, c1, c2; + + /* Suppress warning about unused parameters. */ + (void) np; (void) pmobf; (void) gravcap_f; (void) fflux; + + /* + * Define face pressures as weighted average of connecting cell + * pressures. Specifically, we define + * + * pf = (t1 p1 + t2 p2) / (t1 + t2) + * + * in which t1 and t2 are the one-sided transmissibilities and p1 + * and p2 are the associated cell pressures. + * + * NOTE: The formula does not account for effects of gravity or + * flux boundary conditions. + */ + for (f = 0; f < G->number_of_faces; f++) { + scratch_f[f] = fpress[f] = 0.0; + } + + for (c = i = 0; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + scratch_f[f] += htrans[i]; + fpress[f] += htrans[i] * cpress[c]; + } + } + + for (f = 0; f < G->number_of_faces; f++) { + fpress[f] /= scratch_f[f]; + + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; + + if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == PRESSURE)) { + fpress[f] = bc->bcval[f]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_flux(grid_t *G, + flowbc_t *bc, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + const double *cpress, + double *fflux) +/* ---------------------------------------------------------------------- */ +{ + int f, c1, c2, p; + double t, dp, g; + + for (f = 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] == FLUX)) { + fflux[f] = bc->bcval[f]; + continue; + } + + t = g = 0.0; + for (p = 0; p < np; p++) { + t += pmobf[f*np + p]; + g += pmobf[f*np + p] * gravcap_f[f*np + p]; + } + /* t *= trans[f]; */ + + if ((c1 >= 0) && (c2 >= 0)) { + dp = cpress[c1] - cpress[c2]; + } else if (bc->type[f] == PRESSURE) { + if (c1 < 0) { + dp = bc->bcval[f] - cpress[c2]; + /* g = -g; */ + } else { + dp = cpress[c1] - bc->bcval[f]; + } + } else { + /* No BC -> no-flow (== pressure drop offsets gravity) */ + dp = -g / t; + } + + fflux[f] = trans[f] * (t*dp + g); + } +} + + +static size_t +maxconn(grid_t *G) +{ + int c; + size_t m, n; + + for (c = 0, m = 0; c < G->number_of_cells; c++) { + n = G->cell_facepos[c + 1] - G->cell_facepos[c]; + + if (n > m) { m = n; } + } + + return m; +} + + +/* ---------------------------------------------------------------------- */ +static int +is_incompr(int nc, const double *accum) +/* ---------------------------------------------------------------------- */ +{ + int c, incompr; + + for (c = 0, incompr = 1; (c < nc) && incompr; c++) { + incompr = ! (fabs(accum[c]) > 0.0); + } + + return incompr; +} + + +/* ====================================================================== + * Public interface below separator. + * ====================================================================== */ + +/* ---------------------------------------------------------------------- */ +struct cfs_tpfa_data * +cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) +/* ---------------------------------------------------------------------- */ +{ + size_t nc, nf, ngconn, nwconn; + struct cfs_tpfa_data *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->pimpl = impl_allocate(G, W, maxconn(G), nphases); + new->A = construct_matrix(G, W); + + if ((new->pimpl == NULL) || (new->A == NULL)) { + cfs_tpfa_destroy(new); + new = NULL; + } + } + + if (new != NULL) { + nc = G->number_of_cells; + nf = G->number_of_faces; + ngconn = G->cell_facepos[nc]; + nwconn = 0; + + if (W != NULL) { + nwconn = W->well_connpos[ W->number_of_wells ]; + } + + /* Allocate linear system components */ + 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); + + new->pimpl->flux_work = + new->pimpl->compflux_deriv_f + (nphases * 2 * nf); + + new->pimpl->scratch_f = + new->pimpl->flux_work + (nphases * (1 + 2)); + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_assemble(grid_t *G, + double dt, + flowbc_t *bc, + const double *src, + const double *zc, + struct compr_quantities *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *porevol, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int res_is_neumann, c, np2; + + csrmatrix_zero( h->A); + vector_zero (h->A->m, h->b); + + h->pimpl->is_incomp = 1; + + 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)) { + assert (wdata != NULL); + well_is_neumann = assemble_well_contrib(G->number_of_cells, + W, wctrl, h); + } 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) {*/ + h->A->sa[0] *= 2; + } +} + + +#if 0 +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_press_increment(grid_t *G, + well_t *W, + struct cfs_tpfa_data *h, + double *cpress_inc, + double *wpress_inc) +/* ---------------------------------------------------------------------- */ +{ + /* Assign cell pressure directly from solution vector */ + memcpy(cpress_inc, h->x, G->number_of_cells * sizeof *cpress_inc); + + if (W != NULL) { + assert (wpress_inc != NULL); + + /* Assign well BHP directly from solution vector */ + memcpy(wpress_inc, h->x + G->number_of_cells, + W->number_of_wells * sizeof *wpress_inc); + } +} + + +/* ---------------------------------------------------------------------- */ +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, + const double *wpress, + struct completion_data *wdata, + double *fflux, + double *wflux) +/* ---------------------------------------------------------------------- */ +{ + compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux); + + if ((W != NULL) && (wdata != NULL)) { + assert (wpress != NULL); + assert (wflux != NULL); + + compute_wflux(W, np, wdata, cpress, wpress, wflux); + } +} + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_data *h, + const double *cpress, + const double *fflux, + double *fpress) +/* ---------------------------------------------------------------------- */ +{ + compute_fpress(G, bc, np, htrans, pmobf, gravcap_f, + cpress, fflux, fpress, h->pimpl->scratch_f); +} + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_retrieve_masstrans(grid_t *G, + int np, + struct cfs_tpfa_data *h, + double *masstrans_f) +/* ---------------------------------------------------------------------- */ +{ + memcpy(masstrans_f, h->pimpl->masstrans_f, + np * G->number_of_faces * sizeof *masstrans_f); +} + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_retrieve_gravtrans(grid_t *G, + int np, + struct cfs_tpfa_data *h, + double *gravtrans_f) +/* ---------------------------------------------------------------------- */ +{ + memcpy(gravtrans_f, h->pimpl->gravtrans_f, + np * G->number_of_faces * sizeof *gravtrans_f); +} + + +/* ---------------------------------------------------------------------- */ +static double +cfs_tpfa_impes_maxtime_cell(int c, + grid_t *G, + struct compr_quantities *cq, + const double *trans, + const double *porevol, + struct cfs_tpfa_data *h, + const double *dpmobf, + const double *surf_dens, + const double *gravity) +/* ---------------------------------------------------------------------- */ +{ + /* Reference: + K. H. Coats, "IMPES Stability: The Stable Step", SPE 69225 + + Capillary pressure parts not included. + */ + + int i, j, k, f, c2; + double f11, f12, f21, f22; + double dp, dzg, tr, tmob, detF, eqv_flux; + const double *pmob; + const double *A; + /* This is intended to be compatible with the dune-porsol blackoil + code. This library is otherwise not depending on particular + orderings of phases or components, so at some point this + function should be generalized. */ + enum { Water = 0, Oil = 1, Gas = 2 }; + enum { num_phases = 3 }; + double rho[num_phases]; + double pot[num_phases]; + /* Notation: dpmob[Oil][Water] is d/ds_w(lambda_o) */ + double dpmob[num_phases][num_phases] + = { {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0} }; + + assert (cq->nphases == 3); + + f11 = f12 = f21 = f22 = 0.0; + + /* Loop over neighbour faces to accumulate f11, f12 etc. */ + for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; ++i) { + f = G->cell_faces[i]; + if ((c2 = G->face_cells[2*f + 0]) == c) { + c2 = G->face_cells[2*f + 1]; + } + + /* Initially only interior faces */ + if (c2 < 0) { + continue; + } + + /* Computing density */ + A = cq->Af + f*(cq->nphases)*(cq->nphases); + for (j = 0; j < cq->nphases; ++j) { + rho[j] = 0.0; + for (k = 0; k < cq->nphases; ++k) { + rho[j] += A[cq->nphases*j + k]*surf_dens[k]; + } + } + /* Computing gravity potentials */ + dp = h->x[c] - h->x[c2]; + dzg = 0.0; + for (j = 0; j < G->dimensions; ++j) { + dzg += (G->cell_centroids[G->dimensions*c + j] - G->cell_centroids[G->dimensions*c2 + j])*gravity[j]; + } + for (j = 0; j < cq->nphases; ++j) { + pot[j] = fabs(dp - rho[j]*dzg); + } + /* Filling the dpmob array from available data. + Note that we only need the following combinations: + (Water, Water) + (Water, Gas) + (Oil, Water) + (Oil, Gas) + (Gas, Gas) + + No derivatives w.r.t. Oil is needed, since there are only two + independent saturation variables. + + The lack of (Gas, Water) may be due to assumptions on the + three-phase model used (should be checked to be compatible + with our choices). + */ + dpmob[Water][Water] = dpmobf[9*f]; + dpmob[Water][Gas] = dpmobf[9*f + 2]; + dpmob[Oil][Water] = dpmobf[9*f + 3]; + dpmob[Oil][Gas] = dpmobf[9*f + 5]; + dpmob[Gas][Gas] = dpmobf[9*f + 8]; + /* Computing the flux parts f_ij */ + pmob = cq->phasemobf + f*cq->nphases; + tr = trans[f]; + tmob = pmob[Water] + pmob[Oil] + pmob[Gas]; + f11 += tr*((pmob[Oil] + pmob[Gas])*dpmob[Water][Water]*pot[Water] + - pmob[Water]*dpmob[Oil][Water]*pot[Oil])/tmob; + f12 += -tr*(pmob[Water]*dpmob[Oil][Gas]*pot[Oil] + + pmob[Water]*dpmob[Gas][Gas]*pot[Gas] + - (pmob[Oil] + pmob[Gas])*dpmob[Water][Gas]*pot[Water])/tmob; + f21 += -tr*(pmob[Gas]*dpmob[Water][Water]*pot[Water] + + pmob[Gas]*dpmob[Oil][Water]*pot[Oil])/tmob; + f22 += tr*(-pmob[Gas]*dpmob[Oil][Gas]*pot[Oil] + + (pmob[Water] + pmob[Oil])*dpmob[Gas][Gas]*pot[Gas] + - pmob[Gas]*dpmob[Water][Gas]*pot[Water])/tmob; + } + + /* (from eq. 3, 4a-e, 5a-c) + F_i = 1/2 |f11_i + f22_i + \sqrt{G}| + G = (f11_i + f22_i)^2 - 4 det(F_i) + fXX_i = \sum_j fXX_ij (j runs over the neighbours) + det(F_i) = f11_i f22_i - f12_i f21_i + */ + detF = f11*f22 - f12*f21; + eqv_flux = 0.5*fabs(f11 + f22 + sqrt((f11 + f22)*(f11 + f22) - 4*detF)); + /* delta_t < porevol/eqv_flux */ + return porevol[c]/eqv_flux; +} + +/* ---------------------------------------------------------------------- */ +double +cfs_tpfa_impes_maxtime(grid_t *G, + struct compr_quantities *cq, + const double *trans, + const double *porevol, + struct cfs_tpfa_data *h, + const double *dpmobf, + const double *surf_dens, + const double *gravity) +/* ---------------------------------------------------------------------- */ +{ + int c; + double max_dt, cell_dt; + max_dt = 1e100; + for (c = 0; c < G->number_of_cells; ++c) { + cell_dt = cfs_tpfa_impes_maxtime_cell(c, G, cq, trans, porevol, h, + dpmobf, surf_dens, gravity); + if (cell_dt < max_dt) { + max_dt = cell_dt; + } + } + return max_dt; +} + + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_expl_mass_transport(grid_t *G, + well_t *W, + struct completion_data *wdata, + int np, + double dt, + const double *porevol, + struct cfs_tpfa_data *h, + double *surf_vol) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f, c2, p, w; + double dp, dz, gsgn; + const double *masstrans_f, *gravtrans_f, *masstrans_p, *gravtrans_p; + const double *cpress, *wpress; + + /* Suppress warning about unused parameter. */ + (void) wdata; + + masstrans_f = h->pimpl->masstrans_f; + gravtrans_f = h->pimpl->gravtrans_f; + masstrans_p = h->pimpl->masstrans_p; + gravtrans_p = h->pimpl->gravtrans_p; + cpress = h->x; + wpress = h->x + G->number_of_cells; + + /* Transport through interior faces */ + for (c = i = 0; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + + if ((c2 = G->face_cells[2*f + 0]) == c) { + gsgn = 1.0; + c2 = G->face_cells[2*f + 1]; + } else { + gsgn = -1.0; + } + + if (c2 >= 0) { + dp = cpress[c] - cpress[c2]; + + for (p = 0; p < np; p++) { + dz = masstrans_f[f*np + p] * dp; + dz += gravtrans_f[f*np + p] * gsgn; + + /* dz > 0 => flow from c into c2. */ + surf_vol[c*np + p] -= dz * dt / porevol[c]; + } + } + } + } + + /* Transport through well perforations */ + if (W != NULL) { + 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]; + + for (p = 0; p < np; p++) { + dz = masstrans_p[i*np + p] * dp; + dz += gravtrans_p[i*np + p]; + + /* dz > 0 => flow from perforation into c. */ + surf_vol[c*np + p] += dz * dt / porevol[c]; + } + } + } + } +} +#endif + + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_destroy(struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + if (h != NULL) { + csrmatrix_delete(h->A); + impl_deallocate (h->pimpl); + } + + free(h); +} diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h new file mode 100644 index 00000000..e96c1e93 --- /dev/null +++ b/src/cfs_tpfa_residual.h @@ -0,0 +1,128 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_CFS_TPFA_HEADER_INCLUDED +#define OPM_CFS_TPFA_HEADER_INCLUDED + +#include "grid.h" +#include "flow_bc.h" +#include "well.h" + +#ifdef __cplusplus +extern "C" { +#endif + +struct cfs_tpfa_impl; +struct CSRMatrix; +struct compr_quantities; + +struct cfs_tpfa_data { + struct CSRMatrix *A; + double *b; + double *x; + + struct cfs_tpfa_impl *pimpl; +}; + + +struct cfs_tpfa_data * +cfs_tpfa_construct(grid_t *G, well_t *W, int nphases); + +void +cfs_tpfa_assemble(grid_t *G, + double dt, + flowbc_t *bc, + const double *src, + const double *zc, + struct compr_quantities *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *porevol, + struct cfs_tpfa_data *h); + +void +cfs_tpfa_press_increment(grid_t *G, + struct cfs_tpfa_data *h, + double *cpress_inc); + +#if 0 +void +cfs_tpfa_flux(grid_t *G, + flowbc_t *bc, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + const double *cpress, + double *fflux); +#endif + +void +cfs_tpfa_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_data *h, + const double *cpress, + const double *fflux, + double *fpress); + +void +cfs_tpfa_retrieve_masstrans(grid_t *G, + int np, + struct cfs_tpfa_data *h, + double *masstrans_f); + +void +cfs_tpfa_retrieve_gravtrans(grid_t *G, + int np, + struct cfs_tpfa_data *h, + double *gravtrans_f); + +double +cfs_tpfa_impes_maxtime(grid_t *G, + struct compr_quantities *cq, + const double *trans, + const double *porevol, + struct cfs_tpfa_data *h, + const double *dpmobf, + const double *surf_dens, + const double *gravity); + +void +cfs_tpfa_expl_mass_transport(grid_t *G, + well_t *W, + struct completion_data *wdata, + int np, + double dt, + const double *porevol, + struct cfs_tpfa_data *h, + double *surf_vol); + +void +cfs_tpfa_destroy(struct cfs_tpfa_data *h); + +#ifdef __cplusplus +} +#endif + +#endif /* OPM_CFS_TPFA_HEADER_INCLUDED */ From 8a63636c94b0b17ce917ee63e13175a114343272 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 08/36] 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 8ce3efc0..0fea19f9 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 e96c1e93..7cbec29a 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, From 986e8183c4b6657a175b49592505bc690c7c0fe1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:17:02 +0200 Subject: [PATCH 09/36] Backed out changeset 6b1a90716ea5 --- src/compr_quant.c | 52 ----------------------------------------------- src/compr_quant.h | 12 +++-------- 2 files changed, 3 insertions(+), 61 deletions(-) diff --git a/src/compr_quant.c b/src/compr_quant.c index 93d12a75..296168c7 100644 --- a/src/compr_quant.c +++ b/src/compr_quant.c @@ -17,61 +17,9 @@ along with OPM. If not, see . */ -#include - -#include "sparse_sys.h" #include "compr_quant.h" -void -compr_quantities_deallocate(struct compr_quantities *cq) -{ - if (cq != NULL) { - free(cq->Ac); - } - - free(cq); -} - - -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np) -{ - size_t alloc_sz, np2; - struct compr_quantities *cq; - - cq = malloc(1 * sizeof *cq); - - if (cq != NULL) { - np2 = np * np; - - alloc_sz = np2 * nc; /* Ac */ - alloc_sz += np2 * nc; /* dAc */ - alloc_sz += np2 * nf; /* Af */ - alloc_sz += np * nf; /* phasemobf */ - alloc_sz += 1 * nc; /* voldiscr */ - - cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); - - if (cq->Ac == NULL) { - compr_quantities_deallocate(cq); - cq = NULL; - } else { - cq->dAc = cq->Ac + (np2 * nc); - cq->Af = cq->dAc + (np2 * nc); - cq->phasemobf = cq->Af + (np2 * nf); - cq->voldiscr = cq->phasemobf + (np * nf); - - cq->nphases = np; - - vector_zero(alloc_sz, cq->Ac); - } - } - - return cq; -} - - /* ---------------------------------------------------------------------- */ /* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ /* ---------------------------------------------------------------------- */ diff --git a/src/compr_quant.h b/src/compr_quant.h index b3531176..cacfb551 100644 --- a/src/compr_quant.h +++ b/src/compr_quant.h @@ -29,21 +29,15 @@ extern "C" { #endif struct compr_quantities { - int nphases; /* Number of phases/components */ + int nphases; /* Number of phases/components */ + double *totcompr; /* Total compressibility per cell */ + double *voldiscr; /* Volume discrepancy per cell */ double *Ac; /* RB^{-1} per cell */ - double *dAc; /* d/dp (RB^{-1}) per cell */ double *Af; /* RB^{-1} per face */ double *phasemobf; /* Phase mobility per face */ - double *voldiscr; /* Volume discrepancy per cell */ }; -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np); - -void -compr_quantities_deallocate(struct compr_quantities *cq); - void compr_flux_term(grid_t *G, const double *fflux, From e54f1cf1475098dfd84862936b759920100966a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:18:22 +0200 Subject: [PATCH 10/36] Backed out changeset 3b7e5d602aab We need a copy of 'compr_quant' before restoring compat. --- src/compr_quant.c | 52 +++++++++++++++++++++++++++++++++++++++++++++++ src/compr_quant.h | 12 ++++++++--- 2 files changed, 61 insertions(+), 3 deletions(-) diff --git a/src/compr_quant.c b/src/compr_quant.c index 296168c7..93d12a75 100644 --- a/src/compr_quant.c +++ b/src/compr_quant.c @@ -17,9 +17,61 @@ along with OPM. If not, see . */ +#include + +#include "sparse_sys.h" #include "compr_quant.h" +void +compr_quantities_deallocate(struct compr_quantities *cq) +{ + if (cq != NULL) { + free(cq->Ac); + } + + free(cq); +} + + +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np) +{ + size_t alloc_sz, np2; + struct compr_quantities *cq; + + cq = malloc(1 * sizeof *cq); + + if (cq != NULL) { + np2 = np * np; + + alloc_sz = np2 * nc; /* Ac */ + alloc_sz += np2 * nc; /* dAc */ + alloc_sz += np2 * nf; /* Af */ + alloc_sz += np * nf; /* phasemobf */ + alloc_sz += 1 * nc; /* voldiscr */ + + cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); + + if (cq->Ac == NULL) { + compr_quantities_deallocate(cq); + cq = NULL; + } else { + cq->dAc = cq->Ac + (np2 * nc); + cq->Af = cq->dAc + (np2 * nc); + cq->phasemobf = cq->Af + (np2 * nf); + cq->voldiscr = cq->phasemobf + (np * nf); + + cq->nphases = np; + + vector_zero(alloc_sz, cq->Ac); + } + } + + return cq; +} + + /* ---------------------------------------------------------------------- */ /* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ /* ---------------------------------------------------------------------- */ diff --git a/src/compr_quant.h b/src/compr_quant.h index cacfb551..b3531176 100644 --- a/src/compr_quant.h +++ b/src/compr_quant.h @@ -29,15 +29,21 @@ extern "C" { #endif struct compr_quantities { - int nphases; /* Number of phases/components */ + int nphases; /* Number of phases/components */ - double *totcompr; /* Total compressibility per cell */ - double *voldiscr; /* Volume discrepancy per cell */ double *Ac; /* RB^{-1} per cell */ + double *dAc; /* d/dp (RB^{-1}) per cell */ double *Af; /* RB^{-1} per face */ double *phasemobf; /* Phase mobility per face */ + double *voldiscr; /* Volume discrepancy per cell */ }; +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np); + +void +compr_quantities_deallocate(struct compr_quantities *cq); + void compr_flux_term(grid_t *G, const double *fflux, From 8260b59942279605b296d5e98a66a731758555a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:19:53 +0200 Subject: [PATCH 11/36] Grab copy of 'compr_quant' module. --- src/compr_quant_general.c | 133 ++++++++++++++++++++++++++++++++++++++ src/compr_quant_general.h | 70 ++++++++++++++++++++ 2 files changed, 203 insertions(+) create mode 100644 src/compr_quant_general.c create mode 100644 src/compr_quant_general.h diff --git a/src/compr_quant_general.c b/src/compr_quant_general.c new file mode 100644 index 00000000..93d12a75 --- /dev/null +++ b/src/compr_quant_general.c @@ -0,0 +1,133 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include "sparse_sys.h" +#include "compr_quant.h" + + +void +compr_quantities_deallocate(struct compr_quantities *cq) +{ + if (cq != NULL) { + free(cq->Ac); + } + + free(cq); +} + + +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np) +{ + size_t alloc_sz, np2; + struct compr_quantities *cq; + + cq = malloc(1 * sizeof *cq); + + if (cq != NULL) { + np2 = np * np; + + alloc_sz = np2 * nc; /* Ac */ + alloc_sz += np2 * nc; /* dAc */ + alloc_sz += np2 * nf; /* Af */ + alloc_sz += np * nf; /* phasemobf */ + alloc_sz += 1 * nc; /* voldiscr */ + + cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); + + if (cq->Ac == NULL) { + compr_quantities_deallocate(cq); + cq = NULL; + } else { + cq->dAc = cq->Ac + (np2 * nc); + cq->Af = cq->dAc + (np2 * nc); + cq->phasemobf = cq->Af + (np2 * nf); + cq->voldiscr = cq->phasemobf + (np * nf); + + cq->nphases = np; + + vector_zero(alloc_sz, cq->Ac); + } + } + + return cq; +} + + +/* ---------------------------------------------------------------------- */ +/* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ +/* ---------------------------------------------------------------------- */ +void +compr_flux_term(grid_t *G, + const double *fflux, + const double *zeta, + double *Biv) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f; + double s; + + for (c = i = 0; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + s = 2.0*(c == G->face_cells[2*f + 0]) - 1.0; + + Biv[i] = zeta[c] * s * fflux[f]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +/* Compute P == ct .* pv ./ dt */ +/* ---------------------------------------------------------------------- */ +void +compr_accum_term(size_t nc, + double dt, + const double *porevol, + const double *totcompr, + double *P) +/* ---------------------------------------------------------------------- */ +{ + size_t c; + + for (c = 0; c < nc; c++) { + P[c] = totcompr[c] * porevol[c] / dt; + } +} + + +/* ---------------------------------------------------------------------- */ +/* Add pressure accumulation term (P*p_0) to cell sources */ +/* ---------------------------------------------------------------------- */ +void +compr_src_add_press_accum(size_t nc, + const double *p0, + const double *P, + double *src) +/* ---------------------------------------------------------------------- */ +{ + size_t c; + + for (c = 0; c < nc; c++) { + src[c] += P[c] * p0[c]; + } +} diff --git a/src/compr_quant_general.h b/src/compr_quant_general.h new file mode 100644 index 00000000..b3531176 --- /dev/null +++ b/src/compr_quant_general.h @@ -0,0 +1,70 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_COMPR_QUANT_HEADER_INCLUDED +#define OPM_COMPR_QUANT_HEADER_INCLUDED + +#include + +#include "grid.h" + +#ifdef __cplusplus +extern "C" { +#endif + +struct compr_quantities { + int nphases; /* Number of phases/components */ + + double *Ac; /* RB^{-1} per cell */ + double *dAc; /* d/dp (RB^{-1}) per cell */ + double *Af; /* RB^{-1} per face */ + double *phasemobf; /* Phase mobility per face */ + double *voldiscr; /* Volume discrepancy per cell */ +}; + +struct compr_quantities * +compr_quantities_allocate(size_t nc, size_t nf, int np); + +void +compr_quantities_deallocate(struct compr_quantities *cq); + +void +compr_flux_term(grid_t *G, + const double *fflux, + const double *zeta, + double *Biv); + +void +compr_accum_term(size_t nc, + double dt, + const double *porevol, + const double *totcompr, + double *P); + +void +compr_src_add_press_accum(size_t nc, + const double *p0, + const double *P, + double *src); + +#ifdef __cplusplus +} +#endif + +#endif /* OPM_COMPR_QUANT_HEADER_INCLUDED */ From e3b52c3a44216d0b13d3bf343932a83d90ac1586 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:20:34 +0200 Subject: [PATCH 12/36] Backed out changeset 9faef9e37070 --- src/compr_quant.c | 52 ----------------------------------------------- src/compr_quant.h | 12 +++-------- 2 files changed, 3 insertions(+), 61 deletions(-) diff --git a/src/compr_quant.c b/src/compr_quant.c index 93d12a75..296168c7 100644 --- a/src/compr_quant.c +++ b/src/compr_quant.c @@ -17,61 +17,9 @@ along with OPM. If not, see . */ -#include - -#include "sparse_sys.h" #include "compr_quant.h" -void -compr_quantities_deallocate(struct compr_quantities *cq) -{ - if (cq != NULL) { - free(cq->Ac); - } - - free(cq); -} - - -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np) -{ - size_t alloc_sz, np2; - struct compr_quantities *cq; - - cq = malloc(1 * sizeof *cq); - - if (cq != NULL) { - np2 = np * np; - - alloc_sz = np2 * nc; /* Ac */ - alloc_sz += np2 * nc; /* dAc */ - alloc_sz += np2 * nf; /* Af */ - alloc_sz += np * nf; /* phasemobf */ - alloc_sz += 1 * nc; /* voldiscr */ - - cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); - - if (cq->Ac == NULL) { - compr_quantities_deallocate(cq); - cq = NULL; - } else { - cq->dAc = cq->Ac + (np2 * nc); - cq->Af = cq->dAc + (np2 * nc); - cq->phasemobf = cq->Af + (np2 * nf); - cq->voldiscr = cq->phasemobf + (np * nf); - - cq->nphases = np; - - vector_zero(alloc_sz, cq->Ac); - } - } - - return cq; -} - - /* ---------------------------------------------------------------------- */ /* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ /* ---------------------------------------------------------------------- */ diff --git a/src/compr_quant.h b/src/compr_quant.h index b3531176..cacfb551 100644 --- a/src/compr_quant.h +++ b/src/compr_quant.h @@ -29,21 +29,15 @@ extern "C" { #endif struct compr_quantities { - int nphases; /* Number of phases/components */ + int nphases; /* Number of phases/components */ + double *totcompr; /* Total compressibility per cell */ + double *voldiscr; /* Volume discrepancy per cell */ double *Ac; /* RB^{-1} per cell */ - double *dAc; /* d/dp (RB^{-1}) per cell */ double *Af; /* RB^{-1} per face */ double *phasemobf; /* Phase mobility per face */ - double *voldiscr; /* Volume discrepancy per cell */ }; -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np); - -void -compr_quantities_deallocate(struct compr_quantities *cq); - void compr_flux_term(grid_t *G, const double *fflux, From b70a7d6832f1ec76726b32998ea5b27b9792b790 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 23:22:26 +0200 Subject: [PATCH 13/36] Backed out changeset d83710dd6839 --- src/cfs_tpfa.c | 50 ++++++++++++++++++-------------------------------- src/cfs_tpfa.h | 32 +++++++++++++------------------- 2 files changed, 31 insertions(+), 51 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 0fea19f9..90d58eeb 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -992,48 +992,34 @@ cfs_tpfa_assemble(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_press_increment(grid_t *G, - well_t *W, - struct cfs_tpfa_data *h, - double *cpress_inc, - double *wpress_inc) +cfs_tpfa_press_flux(grid_t *G, + flowbc_t *bc, + well_t *W, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + struct completion_data *wdata, + struct cfs_tpfa_data *h, + double *cpress, + double *fflux, + double *wpress, + double *wflux) /* ---------------------------------------------------------------------- */ { /* Assign cell pressure directly from solution vector */ - memcpy(cpress_inc, h->x, G->number_of_cells * sizeof *cpress_inc); + memcpy(cpress, h->x, G->number_of_cells * sizeof *cpress); - if (W != NULL) { - assert (wpress_inc != NULL); - - /* Assign well BHP directly from solution vector */ - memcpy(wpress_inc, h->x + G->number_of_cells, - W->number_of_wells * sizeof *wpress_inc); - } -} - - -/* ---------------------------------------------------------------------- */ -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, - const double *wpress, - struct completion_data *wdata, - double *fflux, - double *wflux) -/* ---------------------------------------------------------------------- */ -{ compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux); if ((W != NULL) && (wdata != NULL)) { assert (wpress != NULL); assert (wflux != NULL); + /* Assign well BHP directly from solution vector */ + memcpy(wpress, h->x + G->number_of_cells, + W->number_of_wells * sizeof *wpress); + compute_wflux(W, np, wdata, cpress, wpress, wflux); } } diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index 7cbec29a..acd809ea 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -60,25 +60,19 @@ cfs_tpfa_assemble(grid_t *G, 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 *wpress_inc); - -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, - const double *wpress, - struct completion_data *wdata, - double *fflux, - double *wflux); +cfs_tpfa_press_flux(grid_t *G, + flowbc_t *bc, + well_t *W, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + struct completion_data *wdata, + struct cfs_tpfa_data *h, + double *cpress, + double *fflux, + double *wpress, + double *wflux); void cfs_tpfa_fpress(grid_t *G, From 3393aa0d9eaf23bde58692e7606531a160e3ff4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 13:01:06 +0200 Subject: [PATCH 14/36] Hook residual formulation up to build. --- src/Makefile.am | 4 + src/cfs_tpfa_residual.c | 266 ++++++++++++++------------------------ src/cfs_tpfa_residual.h | 67 +++++----- src/compr_quant_general.c | 75 +---------- src/compr_quant_general.h | 29 +---- 5 files changed, 144 insertions(+), 297 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index 7f025915..e0da61de 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -8,7 +8,9 @@ blas_lapack.h \ coarse_conn.h \ coarse_sys.h \ compr_quant.h \ +compr_quant_general.h \ cfs_tpfa.h \ +cfs_tpfa_residual.h \ dfs.h \ flow_bc.h \ fsh.h \ @@ -33,7 +35,9 @@ cfsh.c \ coarse_conn.c \ coarse_sys.c \ compr_quant.c \ +compr_quant_general.c \ cfs_tpfa.c \ +cfs_tpfa_residual.c \ dfs.c \ flow_bc.c \ fsh.c \ diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 8ce3efc0..5af4d60c 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -8,9 +8,9 @@ #include "flow_bc.h" #include "well.h" -#include "compr_quant.h" +#include "compr_quant_general.h" #include "trans_tpfa.h" -#include "cfs_tpfa.h" +#include "cfs_tpfa_residual.h" #include "sparse_sys.h" @@ -35,7 +35,7 @@ struct densrat_util { }; -struct cfs_tpfa_impl { +struct cfs_tpfa_res_impl { int is_incomp; /* One entry per component per face */ @@ -110,7 +110,7 @@ allocate_densrat(size_t max_conn, int np) /* ---------------------------------------------------------------------- */ static void -impl_deallocate(struct cfs_tpfa_impl *pimpl) +impl_deallocate(struct cfs_tpfa_res_impl *pimpl) /* ---------------------------------------------------------------------- */ { if (pimpl != NULL) { @@ -123,12 +123,12 @@ 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) +static struct cfs_tpfa_res_impl * +impl_allocate(grid_t *G, size_t max_conn, int np) /* ---------------------------------------------------------------------- */ { size_t nnu, ngconn, nwperf; - struct cfs_tpfa_impl *new; + struct cfs_tpfa_res_impl *new; size_t ddata_sz; @@ -136,13 +136,8 @@ impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) ngconn = G->cell_facepos[ G->number_of_cells ]; nwperf = 0; - if (W != NULL) { - nnu += W->number_of_wells; - nwperf = W->well_connpos[ W->number_of_wells ]; - } - /* Linear system */ - ddata_sz = 2 * nnu; /* b, x */ + ddata_sz = nnu; /* b */ /* Reservoir */ ddata_sz += np * G->number_of_faces ; /* compflux_f */ @@ -170,18 +165,15 @@ impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) /* ---------------------------------------------------------------------- */ static struct CSRMatrix * -construct_matrix(grid_t *G, well_t *W) +construct_matrix(grid_t *G) /* ---------------------------------------------------------------------- */ { - int f, c1, c2, w, i, nc, nnu; + int f, c1, c2, i, nc, nnu; size_t nnz; struct CSRMatrix *A; nc = nnu = G->number_of_cells; - if (W != NULL) { - nnu += W->number_of_wells; - } A = csrmatrix_new_count_nnz(nnu); @@ -202,18 +194,6 @@ construct_matrix(grid_t *G, well_t *W) } } - if (W != NULL) { - /* Well <-> cell connections */ - for (w = i = 0; w < W->number_of_wells; w++) { - for (; i < W->well_connpos[w + 1]; i++) { - c1 = W->well_cells[i]; - - A->ia[ 0 + c1 + 1 ] += 1; /* c -> w */ - A->ia[ nc + w + 1 ] += 1; /* w -> c */ - } - } - } - nnz = csrmatrix_new_elms_pushback(A); if (nnz == 0) { csrmatrix_delete(A); @@ -238,18 +218,6 @@ construct_matrix(grid_t *G, well_t *W) } } - if (W != NULL) { - /* Fill well <-> cell connections */ - for (w = i = 0; w < W->number_of_wells; w++) { - for (; i < W->well_connpos[w + 1]; i++) { - c1 = W->well_cells[i]; - - A->ja[ A->ia[ 0 + c1 + 1 ] ++ ] = nc + w; - A->ja[ A->ia[ nc + w + 1 ] ++ ] = c1 ; - } - } - } - assert ((size_t) A->ia[ nnu ] == nnz); /* Enforce sorted connection structure per row */ @@ -358,7 +326,7 @@ compute_compflux_and_deriv(grid_t *G , const double *pmobf , const double *gcapf , const double *Af , - struct cfs_tpfa_impl *pimpl ) + struct cfs_tpfa_res_impl *pimpl ) { int c1, c2, f, np2; double dp; @@ -419,7 +387,7 @@ init_cell_contrib(grid_t *G , double pvol , double dt , const double *z , - struct cfs_tpfa_impl *pimpl) + struct cfs_tpfa_res_impl *pimpl) { int c1, c2, f, i, conn, nconn; double *cflx, *dcflx; @@ -469,7 +437,7 @@ compute_cell_contrib(grid_t *G , const double *z , const double *Ac , const double *dAc , - struct cfs_tpfa_impl *pimpl) + struct cfs_tpfa_res_impl *pimpl) { int c1, c2, f, i, nconn, np2, p; MAT_SIZE_T nrhs; @@ -539,7 +507,7 @@ compute_cell_contrib(grid_t *G , static int assemble_cell_contrib(grid_t *G, int c, - struct cfs_tpfa_data *h) + struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ { int c1, c2, i, f, j1, j2; @@ -547,9 +515,9 @@ assemble_cell_contrib(grid_t *G, is_neumann = 1; - j1 = csrmatrix_elm_index(c, c, h->A); + j1 = csrmatrix_elm_index(c, c, h->J); - h->A->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; + h->J->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; for (; i < G->cell_facepos[c + 1]; i++) { f = G->cell_faces[i]; @@ -560,13 +528,13 @@ assemble_cell_contrib(grid_t *G, c2 = (c1 == c) ? c2 : c1; if (c2 >= 0) { - j2 = csrmatrix_elm_index(c, c2, h->A); + j2 = csrmatrix_elm_index(c, c2, h->J); - h->A->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; + h->J->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; } } - h->b[ c ] = h->pimpl->ratio->residual; + h->F[ c ] = h->pimpl->ratio->residual; return 0; } @@ -713,74 +681,83 @@ is_incompr(int nc, const double *accum) * Public interface below separator. * ====================================================================== */ + /* ---------------------------------------------------------------------- */ -struct cfs_tpfa_data * -cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) +void +cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ { - size_t nc, nf, ngconn, nwconn; - struct cfs_tpfa_data *new; + if (h != NULL) { + csrmatrix_delete(h->J); + impl_deallocate (h->pimpl); + } - new = malloc(1 * sizeof *new); + free(h); +} - if (new != NULL) { - new->pimpl = impl_allocate(G, W, maxconn(G), nphases); - new->A = construct_matrix(G, W); - if ((new->pimpl == NULL) || (new->A == NULL)) { - cfs_tpfa_destroy(new); - new = NULL; +/* ---------------------------------------------------------------------- */ +struct cfs_tpfa_res_data * +cfs_tpfa_res_construct(grid_t *G, int nphases) +/* ---------------------------------------------------------------------- */ +{ + size_t nc, nf, ngconn; + struct cfs_tpfa_res_data *h; + + h = malloc(1 * sizeof *h); + + if (h != NULL) { + h->pimpl = impl_allocate(G, maxconn(G), nphases); + h->J = construct_matrix(G); + + if ((h->pimpl == NULL) || (h->J == NULL)) { + cfs_tpfa_res_destroy(h); + h = NULL; } } - if (new != NULL) { + if (h != NULL) { nc = G->number_of_cells; nf = G->number_of_faces; ngconn = G->cell_facepos[nc]; - nwconn = 0; - - if (W != NULL) { - nwconn = W->well_connpos[ W->number_of_wells ]; - } /* Allocate linear system components */ - new->b = new->pimpl->ddata + 0; - new->x = new->b + new->A->m; + h->F = h->pimpl->ddata + 0; - new->pimpl->compflux_f = new->x + new->A->m; - new->pimpl->compflux_deriv_f = - new->pimpl->compflux_f + (nphases * nf); + h->pimpl->compflux_f = h->F + h->J->m; + h->pimpl->compflux_deriv_f = + h->pimpl->compflux_f + (nphases * nf); - new->pimpl->flux_work = - new->pimpl->compflux_deriv_f + (nphases * 2 * nf); + h->pimpl->flux_work = + h->pimpl->compflux_deriv_f + (nphases * 2 * nf); - new->pimpl->scratch_f = - new->pimpl->flux_work + (nphases * (1 + 2)); + h->pimpl->scratch_f = + h->pimpl->flux_work + (nphases * (1 + 2)); } - return new; + return h; } /* ---------------------------------------------------------------------- */ void -cfs_tpfa_assemble(grid_t *G, - double dt, - flowbc_t *bc, - const double *src, - const double *zc, - struct compr_quantities *cq, - const double *trans, - const double *gravcap_f, - const double *cpress, - const double *porevol, - struct cfs_tpfa_data *h) +cfs_tpfa_res_assemble(grid_t *G, + double dt, + flowbc_t *bc, + const double *src, + const double *zc, + struct compr_quantities_gen *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *porevol, + struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ { int res_is_neumann, c, np2; - csrmatrix_zero( h->A); - vector_zero (h->A->m, h->b); + csrmatrix_zero( h->J); + vector_zero (h->J->m, h->F); h->pimpl->is_incomp = 1; @@ -814,70 +791,37 @@ cfs_tpfa_assemble(grid_t *G, if (res_is_neumann && h->pimpl->is_incomp) { /*&& well_is_neumann && h->pimpl->is_incomp) {*/ - h->A->sa[0] *= 2; + h->J->sa[0] *= 2; } } +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_res_flux(grid_t *G, + flowbc_t *bc, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + const double *cpress, + double *fflux) +/* ---------------------------------------------------------------------- */ +{ + compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux); +} + + #if 0 /* ---------------------------------------------------------------------- */ void -cfs_tpfa_press_increment(grid_t *G, - well_t *W, - struct cfs_tpfa_data *h, - double *cpress_inc, - double *wpress_inc) -/* ---------------------------------------------------------------------- */ -{ - /* Assign cell pressure directly from solution vector */ - memcpy(cpress_inc, h->x, G->number_of_cells * sizeof *cpress_inc); - - if (W != NULL) { - assert (wpress_inc != NULL); - - /* Assign well BHP directly from solution vector */ - memcpy(wpress_inc, h->x + G->number_of_cells, - W->number_of_wells * sizeof *wpress_inc); - } -} - - -/* ---------------------------------------------------------------------- */ -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, - const double *wpress, - struct completion_data *wdata, - double *fflux, - double *wflux) -/* ---------------------------------------------------------------------- */ -{ - compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux); - - if ((W != NULL) && (wdata != NULL)) { - assert (wpress != NULL); - assert (wflux != NULL); - - compute_wflux(W, np, wdata, cpress, wpress, wflux); - } -} - - -/* ---------------------------------------------------------------------- */ -void -cfs_tpfa_fpress(grid_t *G, +cfs_tpfa_res_fpress(grid_t *G, flowbc_t *bc, int np, const double *htrans, const double *pmobf, const double *gravcap_f, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, const double *cpress, const double *fflux, double *fpress) @@ -890,9 +834,9 @@ cfs_tpfa_fpress(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_retrieve_masstrans(grid_t *G, +cfs_tpfa_res_retrieve_masstrans(grid_t *G, int np, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, double *masstrans_f) /* ---------------------------------------------------------------------- */ { @@ -903,9 +847,9 @@ cfs_tpfa_retrieve_masstrans(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_retrieve_gravtrans(grid_t *G, +cfs_tpfa_res_retrieve_gravtrans(grid_t *G, int np, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, double *gravtrans_f) /* ---------------------------------------------------------------------- */ { @@ -916,12 +860,12 @@ cfs_tpfa_retrieve_gravtrans(grid_t *G, /* ---------------------------------------------------------------------- */ static double -cfs_tpfa_impes_maxtime_cell(int c, +cfs_tpfa_res_impes_maxtime_cell(int c, grid_t *G, struct compr_quantities *cq, const double *trans, const double *porevol, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, const double *dpmobf, const double *surf_dens, const double *gravity) @@ -1033,11 +977,11 @@ cfs_tpfa_impes_maxtime_cell(int c, /* ---------------------------------------------------------------------- */ double -cfs_tpfa_impes_maxtime(grid_t *G, +cfs_tpfa_res_impes_maxtime(grid_t *G, struct compr_quantities *cq, const double *trans, const double *porevol, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, const double *dpmobf, const double *surf_dens, const double *gravity) @@ -1047,7 +991,7 @@ cfs_tpfa_impes_maxtime(grid_t *G, double max_dt, cell_dt; max_dt = 1e100; for (c = 0; c < G->number_of_cells; ++c) { - cell_dt = cfs_tpfa_impes_maxtime_cell(c, G, cq, trans, porevol, h, + cell_dt = cfs_tpfa_res_impes_maxtime_cell(c, G, cq, trans, porevol, h, dpmobf, surf_dens, gravity); if (cell_dt < max_dt) { max_dt = cell_dt; @@ -1060,13 +1004,13 @@ cfs_tpfa_impes_maxtime(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_expl_mass_transport(grid_t *G, +cfs_tpfa_res_expl_mass_transport(grid_t *G, well_t *W, struct completion_data *wdata, int np, double dt, const double *porevol, - struct cfs_tpfa_data *h, + struct cfs_tpfa_res_data *h, double *surf_vol) /* ---------------------------------------------------------------------- */ { @@ -1130,17 +1074,3 @@ cfs_tpfa_expl_mass_transport(grid_t *G, } } #endif - - -/* ---------------------------------------------------------------------- */ -void -cfs_tpfa_destroy(struct cfs_tpfa_data *h) -/* ---------------------------------------------------------------------- */ -{ - if (h != NULL) { - csrmatrix_delete(h->A); - impl_deallocate (h->pimpl); - } - - free(h); -} diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h index e96c1e93..1a2ed48a 100644 --- a/src/cfs_tpfa_residual.h +++ b/src/cfs_tpfa_residual.h @@ -28,52 +28,48 @@ extern "C" { #endif -struct cfs_tpfa_impl; +struct cfs_tpfa_res_impl; struct CSRMatrix; -struct compr_quantities; +struct compr_quantities_gen; -struct cfs_tpfa_data { - struct CSRMatrix *A; - double *b; - double *x; +struct cfs_tpfa_res_data { + struct CSRMatrix *J; + double *F; - struct cfs_tpfa_impl *pimpl; + struct cfs_tpfa_res_impl *pimpl; }; -struct cfs_tpfa_data * -cfs_tpfa_construct(grid_t *G, well_t *W, int nphases); +struct cfs_tpfa_res_data * +cfs_tpfa_res_construct(grid_t *G, int nphases); void -cfs_tpfa_assemble(grid_t *G, - double dt, - flowbc_t *bc, - const double *src, - const double *zc, - struct compr_quantities *cq, - const double *trans, - const double *gravcap_f, - const double *cpress, - const double *porevol, - struct cfs_tpfa_data *h); +cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); void -cfs_tpfa_press_increment(grid_t *G, - struct cfs_tpfa_data *h, - double *cpress_inc); +cfs_tpfa_res_assemble(grid_t *G, + double dt, + flowbc_t *bc, + const double *src, + const double *zc, + struct compr_quantities_gen *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *porevol, + struct cfs_tpfa_res_data *h); + +void +cfs_tpfa_res_flux(grid_t *G, + flowbc_t *bc, + int np, + const double *trans, + const double *pmobf, + const double *gravcap_f, + const double *cpress, + double *fflux); #if 0 -void -cfs_tpfa_flux(grid_t *G, - flowbc_t *bc, - int np, - const double *trans, - const double *pmobf, - const double *gravcap_f, - const double *cpress, - double *fflux); -#endif - void cfs_tpfa_fpress(grid_t *G, flowbc_t *bc, @@ -118,8 +114,7 @@ cfs_tpfa_expl_mass_transport(grid_t *G, struct cfs_tpfa_data *h, double *surf_vol); -void -cfs_tpfa_destroy(struct cfs_tpfa_data *h); +#endif #ifdef __cplusplus } diff --git a/src/compr_quant_general.c b/src/compr_quant_general.c index 93d12a75..361c1fba 100644 --- a/src/compr_quant_general.c +++ b/src/compr_quant_general.c @@ -20,11 +20,11 @@ #include #include "sparse_sys.h" -#include "compr_quant.h" +#include "compr_quant_general.h" void -compr_quantities_deallocate(struct compr_quantities *cq) +compr_quantities_gen_deallocate(struct compr_quantities_gen *cq) { if (cq != NULL) { free(cq->Ac); @@ -34,11 +34,11 @@ compr_quantities_deallocate(struct compr_quantities *cq) } -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np) +struct compr_quantities_gen * +compr_quantities_gen_allocate(size_t nc, size_t nf, int np) { - size_t alloc_sz, np2; - struct compr_quantities *cq; + size_t alloc_sz, np2; + struct compr_quantities_gen *cq; cq = malloc(1 * sizeof *cq); @@ -54,7 +54,7 @@ compr_quantities_allocate(size_t nc, size_t nf, int np) cq->Ac = malloc(alloc_sz * sizeof *cq->Ac); if (cq->Ac == NULL) { - compr_quantities_deallocate(cq); + compr_quantities_gen_deallocate(cq); cq = NULL; } else { cq->dAc = cq->Ac + (np2 * nc); @@ -70,64 +70,3 @@ compr_quantities_allocate(size_t nc, size_t nf, int np) return cq; } - - -/* ---------------------------------------------------------------------- */ -/* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */ -/* ---------------------------------------------------------------------- */ -void -compr_flux_term(grid_t *G, - const double *fflux, - const double *zeta, - double *Biv) -/* ---------------------------------------------------------------------- */ -{ - int c, i, f; - double s; - - for (c = i = 0; c < G->number_of_cells; c++) { - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; - s = 2.0*(c == G->face_cells[2*f + 0]) - 1.0; - - Biv[i] = zeta[c] * s * fflux[f]; - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Compute P == ct .* pv ./ dt */ -/* ---------------------------------------------------------------------- */ -void -compr_accum_term(size_t nc, - double dt, - const double *porevol, - const double *totcompr, - double *P) -/* ---------------------------------------------------------------------- */ -{ - size_t c; - - for (c = 0; c < nc; c++) { - P[c] = totcompr[c] * porevol[c] / dt; - } -} - - -/* ---------------------------------------------------------------------- */ -/* Add pressure accumulation term (P*p_0) to cell sources */ -/* ---------------------------------------------------------------------- */ -void -compr_src_add_press_accum(size_t nc, - const double *p0, - const double *P, - double *src) -/* ---------------------------------------------------------------------- */ -{ - size_t c; - - for (c = 0; c < nc; c++) { - src[c] += P[c] * p0[c]; - } -} diff --git a/src/compr_quant_general.h b/src/compr_quant_general.h index b3531176..9fd13944 100644 --- a/src/compr_quant_general.h +++ b/src/compr_quant_general.h @@ -22,13 +22,11 @@ #include -#include "grid.h" - #ifdef __cplusplus extern "C" { #endif -struct compr_quantities { +struct compr_quantities_gen { int nphases; /* Number of phases/components */ double *Ac; /* RB^{-1} per cell */ @@ -38,30 +36,11 @@ struct compr_quantities { double *voldiscr; /* Volume discrepancy per cell */ }; -struct compr_quantities * -compr_quantities_allocate(size_t nc, size_t nf, int np); +struct compr_quantities_gen * +compr_quantities_gen_allocate(size_t nc, size_t nf, int np); void -compr_quantities_deallocate(struct compr_quantities *cq); - -void -compr_flux_term(grid_t *G, - const double *fflux, - const double *zeta, - double *Biv); - -void -compr_accum_term(size_t nc, - double dt, - const double *porevol, - const double *totcompr, - double *P); - -void -compr_src_add_press_accum(size_t nc, - const double *p0, - const double *P, - double *src); +compr_quantities_gen_deallocate(struct compr_quantities_gen *cq); #ifdef __cplusplus } From b0645e7c9a3ba37effff4e016ee95af90e08c925 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 13:11:44 +0200 Subject: [PATCH 15/36] Remove further traces of well support. Wells will be introduced at a later time. --- src/cfs_tpfa_residual.c | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 5af4d60c..5331b1ee 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -775,22 +775,9 @@ cfs_tpfa_res_assemble(grid_t *G, assemble_cell_contrib(G, c, h); } -#if 0 - res_is_neumann = assemble_cell_contrib(G, bc, src, h); - - if ((W != NULL) && (wctrl != NULL)) { - assert (wdata != NULL); - well_is_neumann = assemble_well_contrib(G->number_of_cells, - W, wctrl, h); - } 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) {*/ h->J->sa[0] *= 2; } } From cabd8ace6ad5aa2e5e69c8e2ffe203d257f051ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 13:27:24 +0200 Subject: [PATCH 16/36] Re-install interface pressure calculation in residual formulation. --- src/cfs_tpfa_residual.c | 26 +++++++++++++------------- src/cfs_tpfa_residual.h | 22 +++++++++++----------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 5331b1ee..0bbb58cc 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -9,10 +9,10 @@ #include "well.h" #include "compr_quant_general.h" -#include "trans_tpfa.h" -#include "cfs_tpfa_residual.h" #include "sparse_sys.h" +#include "trans_tpfa.h" +#include "cfs_tpfa_residual.h" #if defined(MAX) #undef MAX @@ -799,19 +799,18 @@ cfs_tpfa_res_flux(grid_t *G, } -#if 0 /* ---------------------------------------------------------------------- */ void -cfs_tpfa_res_fpress(grid_t *G, - flowbc_t *bc, - int np, - const double *htrans, - const double *pmobf, - const double *gravcap_f, - struct cfs_tpfa_res_data *h, - const double *cpress, - const double *fflux, - double *fpress) +cfs_tpfa_res_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_res_data *h, + const double *cpress, + const double *fflux, + double *fpress) /* ---------------------------------------------------------------------- */ { compute_fpress(G, bc, np, htrans, pmobf, gravcap_f, @@ -819,6 +818,7 @@ cfs_tpfa_res_fpress(grid_t *G, } +#if 0 /* ---------------------------------------------------------------------- */ void cfs_tpfa_res_retrieve_masstrans(grid_t *G, diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h index 1a2ed48a..afdc639f 100644 --- a/src/cfs_tpfa_residual.h +++ b/src/cfs_tpfa_residual.h @@ -69,19 +69,19 @@ cfs_tpfa_res_flux(grid_t *G, const double *cpress, double *fflux); -#if 0 void -cfs_tpfa_fpress(grid_t *G, - flowbc_t *bc, - int np, - const double *htrans, - const double *pmobf, - const double *gravcap_f, - struct cfs_tpfa_data *h, - const double *cpress, - const double *fflux, - double *fpress); +cfs_tpfa_res_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_res_data *h, + const double *cpress, + const double *fflux, + double *fpress); +#if 0 void cfs_tpfa_retrieve_masstrans(grid_t *G, int np, From 21701c8eff14ca78d8ad96d0d6508ea65884376b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 15:11:21 +0200 Subject: [PATCH 17/36] Correct various indexing errors. --- src/cfs_tpfa_residual.c | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 0bbb58cc..425f6cd3 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -85,7 +85,7 @@ allocate_densrat(size_t max_conn, int np) alloc_sz = np * np; /* lu */ alloc_sz += 2 * np; /* t1, t2 */ - alloc_sz += max_conn * 1 ; /* mat_row */ + alloc_sz += (max_conn + 1) * 1 ; /* mat_row */ alloc_sz += (max_conn + 1) * 1 ; /* coeff */ alloc_sz += n_buffer_col * np; /* linsolve_buffer */ @@ -99,7 +99,7 @@ allocate_densrat(size_t max_conn, int np) 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->coeff = ratio->mat_row + ((max_conn + 1) * 1 ); ratio->linsolve_buffer = ratio->coeff + ((max_conn + 1) * 1 ); } } @@ -402,7 +402,7 @@ init_cell_contrib(grid_t *G , cflx = pimpl->ratio->linsolve_buffer + (1 * np); dcflx = cflx + (nconn * np); - for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; c++) { + 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]; @@ -439,7 +439,7 @@ compute_cell_contrib(grid_t *G , const double *dAc , struct cfs_tpfa_res_impl *pimpl) { - int c1, c2, f, i, nconn, np2, p; + int c1, c2, f, i, off, nconn, np2, p; MAT_SIZE_T nrhs; double s, dF1, dF2, *dv, *dv1, *dv2; @@ -479,8 +479,9 @@ compute_cell_contrib(grid_t *G , 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++) { + dv = pimpl->ratio->linsolve_buffer + (1 + nconn)*np; + off = 0; + for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; i++, off++) { f = G->cell_faces[i]; c1 = G->face_cells[2*f + 0]; @@ -496,8 +497,8 @@ compute_cell_contrib(grid_t *G , dF2 += dv2[ p ]; } - pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; - pimpl->ratio->mat_row[ i ] += s * dt * dF2; + pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; + pimpl->ratio->mat_row[ off + 1 ] += s * dt * dF2; } } } @@ -510,7 +511,7 @@ assemble_cell_contrib(grid_t *G, struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ { - int c1, c2, i, f, j1, j2; + int c1, c2, i, f, j1, j2, off; int is_neumann; is_neumann = 1; @@ -519,7 +520,8 @@ assemble_cell_contrib(grid_t *G, h->J->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; - for (; i < G->cell_facepos[c + 1]; i++) { + off = 1; + for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; i++, off++) { f = G->cell_faces[i]; c1 = G->face_cells[2*f + 0]; @@ -530,7 +532,7 @@ assemble_cell_contrib(grid_t *G, if (c2 >= 0) { j2 = csrmatrix_elm_index(c, c2, h->J); - h->J->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; + h->J->sa[j2] += h->pimpl->ratio->mat_row[ off ]; } } @@ -761,8 +763,8 @@ cfs_tpfa_res_assemble(grid_t *G, h->pimpl->is_incomp = 1; - compute_compflux_and_deriv(G, cq->nphases, cpress, trans, gravcap_f, - cq->phasemobf, cq->Af, h->pimpl); + compute_compflux_and_deriv(G, cq->nphases, cpress, trans, + cq->phasemobf, gravcap_f, cq->Af, h->pimpl); np2 = cq->nphases * cq->nphases; for (c = 0; c < G->number_of_cells; From f14b54cee8a41338d1ea9898454c4a91d4f2131e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 15:12:16 +0200 Subject: [PATCH 18/36] Simplify offset calculation. --- src/cfs_tpfa_residual.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 425f6cd3..cde45f28 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -480,7 +480,7 @@ compute_cell_contrib(grid_t *G , /* Accumulate inter-cell Jacobian contributions */ dv = pimpl->ratio->linsolve_buffer + (1 + nconn)*np; - off = 0; + off = 1; for (i = G->cell_facepos[c]; i < G->cell_facepos[c + 1]; i++, off++) { f = G->cell_faces[i]; @@ -497,8 +497,8 @@ compute_cell_contrib(grid_t *G , dF2 += dv2[ p ]; } - pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; - pimpl->ratio->mat_row[ off + 1 ] += s * dt * dF2; + pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; + pimpl->ratio->mat_row[ off ] += s * dt * dF2; } } } From 0f0844b45e8c51b65709cceaf4cd69bbc88268ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 15:14:11 +0200 Subject: [PATCH 19/36] Clear Jacobian row before assembling local contributions. --- src/cfs_tpfa_residual.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index cde45f28..bfdfd67c 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -465,6 +465,8 @@ compute_cell_contrib(grid_t *G , /* Jacobian row */ + vector_zero(1 + nconn, pimpl->ratio->mat_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); From f87b9184aee938591fdc0d7f98f353c2c4d24f65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 19:44:32 +0200 Subject: [PATCH 20/36] Add a simple representation of a compressible (volume flux) source term. --- src/Makefile.am | 2 + src/compr_source.c | 169 +++++++++++++++++++++++++++++++++++++++++++++ src/compr_source.h | 76 ++++++++++++++++++++ 3 files changed, 247 insertions(+) create mode 100644 src/compr_source.c create mode 100644 src/compr_source.h diff --git a/src/Makefile.am b/src/Makefile.am index e0da61de..5d25dae8 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -9,6 +9,7 @@ coarse_conn.h \ coarse_sys.h \ compr_quant.h \ compr_quant_general.h \ +compr_source.h \ cfs_tpfa.h \ cfs_tpfa_residual.h \ dfs.h \ @@ -36,6 +37,7 @@ coarse_conn.c \ coarse_sys.c \ compr_quant.c \ compr_quant_general.c \ +compr_source.c \ cfs_tpfa.c \ cfs_tpfa_residual.c \ dfs.c \ diff --git a/src/compr_source.c b/src/compr_source.c new file mode 100644 index 00000000..09f928b9 --- /dev/null +++ b/src/compr_source.c @@ -0,0 +1,169 @@ +/*=========================================================================== +// +// File: compr_source.c +// +// Created: 2011-10-19 19:22:24+0200 +// +// Authors: Ingeborg S. Ligaarden +// Jostein R. Natvig +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include + +#include "compr_source.h" + +static int +expand_source_tables(int alloc, struct compr_src *src) +{ + int *c; + double *v, *s; + + c = realloc(src->cell , alloc * 1 * sizeof *c); + v = realloc(src->flux , alloc * 1 * sizeof *v); + s = realloc(src->saturation, alloc * src->nphases * sizeof *s); + + if ((c == NULL) || (v == NULL) || (s == NULL)) { + + free(s); free(v); free(c); + + alloc = 0; + } else { + src->cell = c; src->cpty = alloc; + src->flux = v; src->saturation = s; + } + + return alloc; + +} + + +/* ====================================================================== + * Public methods below separator. + * ====================================================================== */ + + +/* ---------------------------------------------------------------------- */ +struct compr_src * +compr_src_allocate(int np, int nsrc) +/* ---------------------------------------------------------------------- */ +{ + int status; + struct compr_src *src; + + if (np <= 0) { + src = NULL; + } else { + src = malloc(1 * sizeof *src); + + if (src != NULL) { + src->nsrc = 0 ; + src->cpty = 0 ; + src->nphases = np; + + src->cell = NULL; + src->flux = NULL; + src->saturation = NULL; + + if (nsrc > 0) { + status = expand_source_tables(nsrc, src); + } else { + status = 1; + } + } + + if (status <= 0) { + compr_src_deallocate(src); + src = NULL; + } + } + + return src; +} + + +/* ---------------------------------------------------------------------- */ +void +compr_src_deallocate(struct compr_src *src) +/* ---------------------------------------------------------------------- */ +{ + if (src != NULL) { + free(src->saturation); + free(src->flux ); + free(src->cell ); + } + + free(src); +} + + +/* ---------------------------------------------------------------------- */ +int +append_compr_source_term(int c , + int np , + double v , + const double *sat, + struct compr_src *src) +/* ---------------------------------------------------------------------- */ +{ + int alloc, status; + + if (np != src->nphases) { + return -1; + } + + if (src->nsrc == src->cpty) { + alloc = (src->cpty > 0) ? 2 * src->cpty : 1; + status = expand_source_tables(alloc, src); + } else { + assert (src->nsrc < src->cpty); + status = 1; + } + + if (status > 0) { + src->cell[ src->nsrc ] = c; + src->flux[ src->nsrc ] = v; + + memcpy(src->saturation + (src->nsrc * np), sat, + np * sizeof *src->saturation); + + src->nsrc += 1; + } + + return status > 0; +} + + +/* ---------------------------------------------------------------------- */ +void +clear_compr_source_term(struct compr_src *src) +/* ---------------------------------------------------------------------- */ +{ + src->nsrc = 0; +} diff --git a/src/compr_source.h b/src/compr_source.h new file mode 100644 index 00000000..d88c1b55 --- /dev/null +++ b/src/compr_source.h @@ -0,0 +1,76 @@ +/*=========================================================================== +// +// File: compr_source.h +// +// Created: 2011-10-19 19:14:30+0200 +// +// Authors: Ingeborg S. Ligaarden +// Jostein R. Natvig +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_COMPR_SOURCE_H_HEADER +#define OPM_COMPR_SOURCE_H_HEADER + +#ifdef __cplusplus +extern "C" { +#endif + +struct compr_src { + int nsrc; + int cpty; + + int nphases; + + int *cell; + double *flux; + double *saturation; +}; + + +struct compr_src * +compr_src_allocate(int np, int nsrc); + +void +compr_src_deallocate(struct compr_src *src); + +int +append_compr_source_term(int c , + int np , + double v , + const double *sat, + struct compr_src *src); + +void +clear_compr_source_term(struct compr_src *src); + + +#ifdef __cplusplus +} +#endif + +#endif /* OPM_COMPR_SOURCE_H_HEADER */ From c2482fbf12299e10fc2ed3e1b8da5d119d4bfc00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 19:48:46 +0200 Subject: [PATCH 21/36] Remove unused function. --- src/cfs_tpfa_residual.c | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index bfdfd67c..e401dbaa 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -36,7 +36,7 @@ struct densrat_util { struct cfs_tpfa_res_impl { - int is_incomp; + int is_incomp /* One entry per component per face */ double *compflux_f; /* A_{ij} v_{ij} */ @@ -666,21 +666,6 @@ maxconn(grid_t *G) } -/* ---------------------------------------------------------------------- */ -static int -is_incompr(int nc, const double *accum) -/* ---------------------------------------------------------------------- */ -{ - int c, incompr; - - for (c = 0, incompr = 1; (c < nc) && incompr; c++) { - incompr = ! (fabs(accum[c]) > 0.0); - } - - return incompr; -} - - /* ====================================================================== * Public interface below separator. * ====================================================================== */ From 7d097bafac9e2208588059a86f7f7f2035c33927 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 19:54:17 +0200 Subject: [PATCH 22/36] Accept a set of compressible source terms rather than a double*. While here, restore an essential semicolon lost in cset 11881e6cd650. --- src/cfs_tpfa_residual.c | 5 +++-- src/cfs_tpfa_residual.h | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index e401dbaa..2cc10b1d 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -9,6 +9,7 @@ #include "well.h" #include "compr_quant_general.h" +#include "compr_source.h" #include "sparse_sys.h" #include "trans_tpfa.h" @@ -36,7 +37,7 @@ struct densrat_util { struct cfs_tpfa_res_impl { - int is_incomp + int is_incomp; /* One entry per component per face */ double *compflux_f; /* A_{ij} v_{ij} */ @@ -733,7 +734,7 @@ void cfs_tpfa_res_assemble(grid_t *G, double dt, flowbc_t *bc, - const double *src, + struct compr_src *src, const double *zc, struct compr_quantities_gen *cq, const double *trans, diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h index afdc639f..332c4279 100644 --- a/src/cfs_tpfa_residual.h +++ b/src/cfs_tpfa_residual.h @@ -31,6 +31,7 @@ extern "C" { struct cfs_tpfa_res_impl; struct CSRMatrix; struct compr_quantities_gen; +struct compr_src; struct cfs_tpfa_res_data { struct CSRMatrix *J; @@ -50,7 +51,7 @@ void cfs_tpfa_res_assemble(grid_t *G, double dt, flowbc_t *bc, - const double *src, + struct compr_src *src, const double *zc, struct compr_quantities_gen *cq, const double *trans, From 735e0fca958330b6408e3f8b93b57297ba846352 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 20:09:51 +0200 Subject: [PATCH 23/36] Implement compressible (volume flux) source term assembly. --- src/cfs_tpfa_residual.c | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 2cc10b1d..313b7447 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -507,6 +507,20 @@ compute_cell_contrib(grid_t *G , } +static void +assemble_sources(struct compr_src *src, struct cfs_tpfa_res_data *h) +{ + int i; + + for (i = 0; i < src->nsrc; i++) { + assert (src->cell[i] >= 0 ); + assert (((size_t) src->cell[i]) < h->J->m); + + h->F[ src->cell[ i ] ] += src->flux[ i ]; + } +} + + /* ---------------------------------------------------------------------- */ static int assemble_cell_contrib(grid_t *G, @@ -765,6 +779,11 @@ cfs_tpfa_res_assemble(grid_t *G, assemble_cell_contrib(G, c, h); } + if (src != NULL) { + assert (src->nphases == cq->nphases); + assemble_sources(src, h); + } + res_is_neumann = 1; if (res_is_neumann && h->pimpl->is_incomp) { From 803e3d1286c750c3d079e77c06050cf207c83cb0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 20:32:50 +0200 Subject: [PATCH 24/36] Check allocation status in correct scope. --- src/compr_source.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/compr_source.c b/src/compr_source.c index 09f928b9..04d12d4c 100644 --- a/src/compr_source.c +++ b/src/compr_source.c @@ -96,11 +96,11 @@ compr_src_allocate(int np, int nsrc) } else { status = 1; } - } - if (status <= 0) { - compr_src_deallocate(src); - src = NULL; + if (status <= 0) { + compr_src_deallocate(src); + src = NULL; + } } } From f28f1137beb1c1f98c6594862cf0675baf13aebf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 22:10:22 +0200 Subject: [PATCH 25/36] Account for flow-direction signs when accumulating fluxes. --- src/cfs_tpfa_residual.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 313b7447..402832d2 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -418,7 +418,7 @@ init_cell_contrib(grid_t *G , cflx += 1 * np; dcflx += 2 * np; - pimpl->ratio->coeff[ conn++ ] = dt; + pimpl->ratio->coeff[ conn++ ] = dt * (2*(c1 == c) - 1.0); } } @@ -516,7 +516,7 @@ assemble_sources(struct compr_src *src, struct cfs_tpfa_res_data *h) assert (src->cell[i] >= 0 ); assert (((size_t) src->cell[i]) < h->J->m); - h->F[ src->cell[ i ] ] += src->flux[ i ]; + h->F[ src->cell[ i ] ] -= src->flux[ i ]; } } From c90fc0d85d5dc396048782590010aa0049957316 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 19 Oct 2011 22:28:41 +0200 Subject: [PATCH 26/36] Don't pretend that the time-step is always one. --- src/cfs_tpfa_residual.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 402832d2..ae530dd3 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -508,7 +508,9 @@ compute_cell_contrib(grid_t *G , static void -assemble_sources(struct compr_src *src, struct cfs_tpfa_res_data *h) +assemble_sources(double dt , + struct compr_src *src, + struct cfs_tpfa_res_data *h ) { int i; @@ -516,7 +518,7 @@ assemble_sources(struct compr_src *src, struct cfs_tpfa_res_data *h) assert (src->cell[i] >= 0 ); assert (((size_t) src->cell[i]) < h->J->m); - h->F[ src->cell[ i ] ] -= src->flux[ i ]; + h->F[ src->cell[ i ] ] -= dt * src->flux[ i ]; } } @@ -781,7 +783,7 @@ cfs_tpfa_res_assemble(grid_t *G, if (src != NULL) { assert (src->nphases == cq->nphases); - assemble_sources(src, h); + assemble_sources(dt, src, h); } res_is_neumann = 1; From 2835762542f6e8452759e88536ae69cfdb6d94a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Oct 2011 13:59:15 +0200 Subject: [PATCH 27/36] Fix remaining bugs for incompressible flows driven by sources. In particular: - Zero sufficient portion of the ->mat_row to hold all connections of a cell in addition to the accumulation term. - Don't write the residual into the accumulation term of the Jacobian matrix row. --- src/cfs_tpfa_residual.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index ae530dd3..aff12494 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -466,20 +466,20 @@ compute_cell_contrib(grid_t *G , /* Jacobian row */ - vector_zero(1 + nconn, pimpl->ratio->mat_row); + vector_zero(1 + (G->cell_facepos[c + 1] - G->cell_facepos[c]), + pimpl->ratio->mat_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; + 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; + pimpl->ratio->mat_row[ 0 ] = - dF2; /* Accumulate inter-cell Jacobian contributions */ dv = pimpl->ratio->linsolve_buffer + (1 + nconn)*np; From b8f7cc176575dd4c608813bd2ce090b3312eb7b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 21 Oct 2011 15:30:44 +0200 Subject: [PATCH 28/36] Don't use uninitialised variables. --- src/cfs_tpfa_residual.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index aff12494..721c70f7 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -369,6 +369,8 @@ count_internal_conn(grid_t *G, int c) { int c1, c2, f, i, nconn; + nconn = 0; + 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]; From 057d291b3eac12da0e01326a2f6f77eff0bb8b7b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 24 Oct 2011 18:10:34 +0200 Subject: [PATCH 29/36] Add a facility for representing compressible boundary conditions. --- src/Makefile.am | 2 + src/compr_bc.c | 182 ++++++++++++++++++++++++++++++++++++++++++++++++ src/compr_bc.h | 83 ++++++++++++++++++++++ 3 files changed, 267 insertions(+) create mode 100644 src/compr_bc.c create mode 100644 src/compr_bc.h diff --git a/src/Makefile.am b/src/Makefile.am index 5d25dae8..602bfe73 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -7,6 +7,7 @@ libopmpressure_HEADERS = \ blas_lapack.h \ coarse_conn.h \ coarse_sys.h \ +compr_bc.h \ compr_quant.h \ compr_quant_general.h \ compr_source.h \ @@ -35,6 +36,7 @@ libopmpressure_la_SOURCES = \ cfsh.c \ coarse_conn.c \ coarse_sys.c \ +compr_bc.c \ compr_quant.c \ compr_quant_general.c \ compr_source.c \ diff --git a/src/compr_bc.c b/src/compr_bc.c new file mode 100644 index 00000000..2e443de0 --- /dev/null +++ b/src/compr_bc.c @@ -0,0 +1,182 @@ +/*=========================================================================== +// +// File: compr_bc.c +// +// Created: 2011-10-24 16:07:17+0200 +// +// Authors: Ingeborg S. Ligaarden +// Jostein R. Natvig +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include + +#include "compr_bc.h" + +static int +expand_source_tables(int alloc, struct compr_bc *bc) +{ + enum compr_bc_type *t; + int *f; + double *p, *v, *s; + + t = realloc(bc->type , alloc * 1 * sizeof *t); + f = realloc(bc->face , alloc * 1 * sizeof *f); + p = realloc(bc->press , alloc * 1 * sizeof *p); + v = realloc(bc->flux , alloc * 1 * sizeof *v); + s = realloc(bc->saturation, alloc * bc->nphases * sizeof *s); + + if ((t == NULL) || (f == NULL) || + (p == NULL) || (v == NULL) || (s == NULL)) { + + free(s); free(v); free(p); free(f); free(t); + + alloc = 0; + } else { + bc->type = t; bc->face = f; + bc->press = p; bc->flux = v; + bc->saturation = s; + } + + return alloc; + +} + + +/* ====================================================================== + * Public methods below separator. + * ====================================================================== */ + + +/* ---------------------------------------------------------------------- */ +struct compr_bc * +compr_bc_allocate(int np, int nbc) +/* ---------------------------------------------------------------------- */ +{ + int status; + struct compr_bc *bc; + + if (np <= 0) { + bc = NULL; + } else { + bc = malloc(1 * sizeof *bc); + + if (bc != NULL) { + bc->nbc = 0 ; + bc->cpty = 0 ; + bc->nphases = np; + + bc->type = NULL; + bc->face = NULL; + bc->press = NULL; + bc->flux = NULL; + bc->saturation = NULL; + + if (nbc > 0) { + status = expand_source_tables(nbc, bc); + } else { + status = 1; + } + + if (status <= 0) { + compr_bc_deallocate(bc); + bc = NULL; + } + } + } + + return bc; +} + + +/* ---------------------------------------------------------------------- */ +void +compr_bc_deallocate(struct compr_bc *bc) +/* ---------------------------------------------------------------------- */ +{ + if (bc != NULL) { + free(bc->saturation); + free(bc->flux ); + free(bc->press ); + free(bc->face ); + free(bc->type ); + } + + free(bc); +} + + +/* ---------------------------------------------------------------------- */ +int +compr_bc_append(enum compr_bc_type t , + int f , + int np , + double p , + double v , + const double *sat, + struct compr_bc *bc ) +/* ---------------------------------------------------------------------- */ +{ + int alloc, status; + + if (np != bc->nphases) { + return -1; + } + + if (bc->nbc == bc->cpty) { + alloc = (bc->cpty > 0) ? 2 * bc->cpty : 1; + status = expand_source_tables(alloc, bc); + } else { + assert (bc->nbc < bc->cpty); + status = 1; + } + + if (status > 0) { + bc->type [ bc->nbc ] = t; + bc->face [ bc->nbc ] = f; + bc->press[ bc->nbc ] = p; + bc->flux [ bc->nbc ] = v; + + memcpy(bc->saturation + (bc->nbc * np), sat, + np * sizeof *bc->saturation); + + bc->nbc += 1; + } + + return status > 0; +} + + +/* ---------------------------------------------------------------------- */ +void +compr_bc_clear(struct compr_bc *bc) +/* ---------------------------------------------------------------------- */ +{ + bc->nbc = 0; +} diff --git a/src/compr_bc.h b/src/compr_bc.h new file mode 100644 index 00000000..e50ee4a3 --- /dev/null +++ b/src/compr_bc.h @@ -0,0 +1,83 @@ +/*=========================================================================== +// +// File: compr_bc.h +// +// Created: 2011-10-24 15:58:16+0200 +// +// Authors: Ingeborg S. Ligaarden +// Jostein R. Natvig +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_COMPR_BC_H_HEADER +#define OPM_COMPR_BC_H_HEADER + +#ifdef __cplusplus +extern "C" { +#endif + +enum compr_bc_type { BC_PRESSURE, BC_RESV }; + +struct compr_bc { + int nbc; + int cpty; + + int nphases; + + enum compr_bc_type *type; + int *face; + + double *press; + double *flux; + double *saturation; +}; + + +struct compr_bc * +compr_bc_allocate(int np, int nbc); + +void +compr_bc_deallocate(struct compr_bc *bc); + +int +compr_bc_append(enum compr_bc_type type, + int f , + int np , + double p , + double v , + const double *sat , + struct compr_bc *bc ); + +void +compr_bc_clear(struct compr_bc *bc); + + +#ifdef __cplusplus +} +#endif + +#endif /* OPM_COMPR_BC_H_HEADER */ From 68c51d8ddf427f16b6f987afa71b2aa063bd72d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 25 Oct 2011 19:34:32 +0200 Subject: [PATCH 30/36] Guard against 'bc' being NULL. This is in preparation of introducing compressible boundary conditions as represented by the 'compr_bc' module. --- src/cfs_tpfa_residual.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 721c70f7..5517f81a 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -612,7 +612,8 @@ compute_fpress(grid_t *G, c1 = G->face_cells[2*f + 0]; c2 = G->face_cells[2*f + 1]; - if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == PRESSURE)) { + if (((c1 < 0) || (c2 < 0)) && + (bc != NULL) && (bc->type[f] == PRESSURE)) { fpress[f] = bc->bcval[f]; } } @@ -638,7 +639,8 @@ compute_flux(grid_t *G, c1 = G->face_cells[2*f + 0]; c2 = G->face_cells[2*f + 1]; - if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == FLUX)) { + if (((c1 < 0) || (c2 < 0)) && + (bc != NULL) && (bc->type[f] == FLUX)) { fflux[f] = bc->bcval[f]; continue; } @@ -652,7 +654,7 @@ compute_flux(grid_t *G, if ((c1 >= 0) && (c2 >= 0)) { dp = cpress[c1] - cpress[c2]; - } else if (bc->type[f] == PRESSURE) { + } else if ((bc != NULL) && (bc->type[f] == PRESSURE)) { if (c1 < 0) { dp = bc->bcval[f] - cpress[c2]; /* g = -g; */ From 38e0cb08ebfcf4dddcbec76852c467d85cfe7a63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 26 Oct 2011 16:19:34 +0200 Subject: [PATCH 31/36] Remove a comment that alludes to a non-portable technique. --- src/coarse_sys.c | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/coarse_sys.c b/src/coarse_sys.c index b2137d37..17503571 100644 --- a/src/coarse_sys.c +++ b/src/coarse_sys.c @@ -975,13 +975,6 @@ unenumerate_local_dofs(size_t cf, { int *b, *c, i; - /* Note that memset()-ing all of m->loc_fno (g->number_of_faces - * entries) may or may not be faster than this loop. The loop is - * entirely random access to m->loc_fno which ruins locality, but - * typically visits only a small subset of the actual elements. - * - * Profiling/measurement needed. */ - for (b = ct->neighbours + 2*(cf + 0); b != ct->neighbours + 2*(cf + 1); b++) { if (*b >= 0) { From 099fa2b61208ebe1c14ee9c98b2370916974bd90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Nov 2011 11:23:30 +0100 Subject: [PATCH 32/36] Name the well_t and well_control_t structures. --- src/well.h | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/well.h b/src/well.h index fe1c64ae..279f57aa 100644 --- a/src/well.h +++ b/src/well.h @@ -28,17 +28,17 @@ extern "C" { enum well_type { INJECTOR, PRODUCER }; enum well_control { BHP , RATE }; -typedef struct { +struct WellCompletions { int number_of_wells; int *well_connpos; int *well_cells; -} well_t; +}; -typedef struct { +struct WellControls { enum well_type *type; enum well_control *ctrl; double *target; -} well_control_t; +}; struct completion_data { double *WI; /* Productivity index */ @@ -47,6 +47,9 @@ struct completion_data { double *phasemob; /* Phase mobility */ }; +typedef struct WellCompletions well_t; +typedef struct WellControls well_control_t; + int allocate_cell_wells(int nc, well_t *W, int **cwpos, int **cwells); From aff142bd54a91669e486155a325afe71dc27dfee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Nov 2011 13:20:37 +0100 Subject: [PATCH 33/36] Take initial steps towards including wells in assembly. Specifically, accept (and currently ignore), a WellCompletions structure into the constructor, and aggregate all driving forces (source terms, boundary conditions and all well-related structures) into an assembler-specific "force" structure. Accept a pointer to such a structure into the assemble() function. Currently ignored except for source terms. --- src/cfs_tpfa_residual.c | 25 +++++++++++++++++-------- src/cfs_tpfa_residual.h | 25 ++++++++++++++++++++----- 2 files changed, 37 insertions(+), 13 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 5517f81a..d6b658e7 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -125,7 +125,10 @@ impl_deallocate(struct cfs_tpfa_res_impl *pimpl) /* ---------------------------------------------------------------------- */ static struct cfs_tpfa_res_impl * -impl_allocate(grid_t *G, size_t max_conn, int np) +impl_allocate(grid_t *G , + struct WellCompletions *wconn , + size_t max_conn, + int np ) /* ---------------------------------------------------------------------- */ { size_t nnu, ngconn, nwperf; @@ -137,6 +140,11 @@ impl_allocate(grid_t *G, size_t max_conn, int np) ngconn = G->cell_facepos[ G->number_of_cells ]; nwperf = 0; + if (wconn != NULL) { + nnu += wconn->number_of_wells; + nwperf = wconn->well_connpos[ wconn->number_of_wells ]; + } + /* Linear system */ ddata_sz = nnu; /* b */ @@ -708,7 +716,9 @@ cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ struct cfs_tpfa_res_data * -cfs_tpfa_res_construct(grid_t *G, int nphases) +cfs_tpfa_res_construct(grid_t *G , + struct WellCompletions *wconn , + int nphases) /* ---------------------------------------------------------------------- */ { size_t nc, nf, ngconn; @@ -717,7 +727,7 @@ cfs_tpfa_res_construct(grid_t *G, int nphases) h = malloc(1 * sizeof *h); if (h != NULL) { - h->pimpl = impl_allocate(G, maxconn(G), nphases); + h->pimpl = impl_allocate(G, wconn, maxconn(G), nphases); h->J = construct_matrix(G); if ((h->pimpl == NULL) || (h->J == NULL)) { @@ -753,8 +763,7 @@ cfs_tpfa_res_construct(grid_t *G, int nphases) void cfs_tpfa_res_assemble(grid_t *G, double dt, - flowbc_t *bc, - struct compr_src *src, + struct cfs_tpfa_res_forces *forces, const double *zc, struct compr_quantities_gen *cq, const double *trans, @@ -785,9 +794,9 @@ cfs_tpfa_res_assemble(grid_t *G, assemble_cell_contrib(G, c, h); } - if (src != NULL) { - assert (src->nphases == cq->nphases); - assemble_sources(dt, src, h); + if ((forces != NULL) && (forces->src != NULL)) { + assert (forces->src->nphases == cq->nphases); + assemble_sources(dt, forces->src, h); } res_is_neumann = 1; diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h index 332c4279..ac75e455 100644 --- a/src/cfs_tpfa_residual.h +++ b/src/cfs_tpfa_residual.h @@ -21,8 +21,6 @@ #define OPM_CFS_TPFA_HEADER_INCLUDED #include "grid.h" -#include "flow_bc.h" -#include "well.h" #ifdef __cplusplus extern "C" { @@ -32,6 +30,22 @@ struct cfs_tpfa_res_impl; struct CSRMatrix; struct compr_quantities_gen; struct compr_src; +struct compr_bc; +struct WellCompletions; +struct WellControls; +struct completion_data; + +struct cfs_tpfa_res_wells { + struct WellCompletions *conn; + struct WellControls *ctrl; + struct completion_data *data; +}; + +struct cfs_tpfa_res_forces { + struct cfs_tpfa_res_wells *W ; + struct compr_bc *bc ; + struct compr_src *src; +}; struct cfs_tpfa_res_data { struct CSRMatrix *J; @@ -42,7 +56,9 @@ struct cfs_tpfa_res_data { struct cfs_tpfa_res_data * -cfs_tpfa_res_construct(grid_t *G, int nphases); +cfs_tpfa_res_construct(grid_t *G , + struct WellCompletions *wconn , + int nphases); void cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); @@ -50,8 +66,7 @@ cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); void cfs_tpfa_res_assemble(grid_t *G, double dt, - flowbc_t *bc, - struct compr_src *src, + struct cfs_tpfa_res_forces *forces, const double *zc, struct compr_quantities_gen *cq, const double *trans, From 8937b6b283836f6074fdadc8e8cede0f026debd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Nov 2011 15:25:48 +0100 Subject: [PATCH 34/36] Allocate backing store for well completion data. Also, include well connections in system matrix. --- src/cfs_tpfa_residual.c | 57 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 51 insertions(+), 6 deletions(-) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index d6b658e7..698dc496 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -43,6 +43,10 @@ struct cfs_tpfa_res_impl { double *compflux_f; /* A_{ij} v_{ij} */ double *compflux_deriv_f; /* A_{ij} \partial_{p} v_{ij} */ + /* One entry per component per perforation */ + double *compflux_p; /* A_{wi} q_{wi} */ + double *compflux_deriv_p; /* A_{wi} \partial_{p} q_{wi} */ + double *flux_work; /* Scratch array for face pressure calculation */ @@ -152,6 +156,10 @@ impl_allocate(grid_t *G , ddata_sz += np * G->number_of_faces ; /* compflux_f */ ddata_sz += np * (2 * G->number_of_faces); /* compflux_deriv_f */ + /* Well perforations */ + ddata_sz += np * nwperf ; /* compflux_p */ + ddata_sz += np * (2 * nwperf); /* compflux_deriv_p */ + ddata_sz += np * (1 + 2) ; /* flux_work */ ddata_sz += 1 * G->number_of_faces ; /* scratch_f */ @@ -174,15 +182,18 @@ impl_allocate(grid_t *G , /* ---------------------------------------------------------------------- */ static struct CSRMatrix * -construct_matrix(grid_t *G) +construct_matrix(grid_t *G, struct WellCompletions *wconn) /* ---------------------------------------------------------------------- */ { - int f, c1, c2, i, nc, nnu; + int f, c1, c2, w, i, nc, nnu; size_t nnz; struct CSRMatrix *A; nc = nnu = G->number_of_cells; + if (wconn != NULL) { + nnu += wconn->number_of_wells; + } A = csrmatrix_new_count_nnz(nnu); @@ -203,6 +214,18 @@ construct_matrix(grid_t *G) } } + if (wconn != NULL) { + /* Well <-> cell connections */ + for (w = i = 0; w < wconn->number_of_wells; w++) { + for (; i < wconn->well_connpos[w + 1]; i++) { + c1 = wconn->well_cells[i]; + + A->ia[ 0 + c1 + 1 ] += 1; /* c -> w */ + A->ia[ nc + w + 1 ] += 1; /* w -> c */ + } + } + } + nnz = csrmatrix_new_elms_pushback(A); if (nnz == 0) { csrmatrix_delete(A); @@ -227,6 +250,18 @@ construct_matrix(grid_t *G) } } + if (wconn != NULL) { + /* Fill well <-> cell connections */ + for (w = i = 0; w < wconn->number_of_wells; w++) { + for (; i < wconn->well_connpos[w + 1]; i++) { + c1 = wconn->well_cells[i]; + + A->ja[ A->ia[ 0 + c1 + 1 ] ++ ] = nc + w; + A->ja[ A->ia[ nc + w + 1 ] ++ ] = c1 ; + } + } + } + assert ((size_t) A->ia[ nnu ] == nnz); /* Enforce sorted connection structure per row */ @@ -721,14 +756,14 @@ cfs_tpfa_res_construct(grid_t *G , int nphases) /* ---------------------------------------------------------------------- */ { - size_t nc, nf, ngconn; + size_t nc, nf, nwperf, ngconn; struct cfs_tpfa_res_data *h; h = malloc(1 * sizeof *h); if (h != NULL) { h->pimpl = impl_allocate(G, wconn, maxconn(G), nphases); - h->J = construct_matrix(G); + h->J = construct_matrix(G, wconn); if ((h->pimpl == NULL) || (h->J == NULL)) { cfs_tpfa_res_destroy(h); @@ -739,18 +774,28 @@ cfs_tpfa_res_construct(grid_t *G , if (h != NULL) { nc = G->number_of_cells; nf = G->number_of_faces; + nwperf = 0; ngconn = G->cell_facepos[nc]; + if (wconn != NULL) { + nwperf = wconn->well_connpos[ wconn->number_of_wells ]; + } + /* Allocate linear system components */ h->F = h->pimpl->ddata + 0; h->pimpl->compflux_f = h->F + h->J->m; - h->pimpl->compflux_deriv_f = + h->pimpl->compflux_p = h->pimpl->compflux_f + (nphases * nf); - h->pimpl->flux_work = + h->pimpl->compflux_deriv_f = + h->pimpl->compflux_p + (nphases * nwperf); + h->pimpl->compflux_deriv_p = h->pimpl->compflux_deriv_f + (nphases * 2 * nf); + h->pimpl->flux_work = + h->pimpl->compflux_deriv_p + (nphases * 2 * nwperf); + h->pimpl->scratch_f = h->pimpl->flux_work + (nphases * (1 + 2)); } From 8677bc3827155dcc83a67d68618cd87a895c0a29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Nov 2011 19:25:47 +0100 Subject: [PATCH 35/36] Compute component fluxes across completions if the model includes wells. This necessitates a public interface change whereby the caller is required to pass the current well (bottom-hole) pressure values into the assembler. --- src/cfs_tpfa_residual.c | 58 +++++++++++++++++++++++++++++++++++++++++ src/cfs_tpfa_residual.h | 1 + 2 files changed, 59 insertions(+) diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c index 698dc496..23d70cfd 100644 --- a/src/cfs_tpfa_residual.c +++ b/src/cfs_tpfa_residual.c @@ -407,6 +407,58 @@ compute_compflux_and_deriv(grid_t *G , } +static void +compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *W , + int np , + const double *cpress, + const double *wpress, + struct cfs_tpfa_res_impl *pimpl ) +{ + int c, i, w, np2; + double pw, dp; + double *WI, *gpot, *Ap, *pmobp; + + double *pflux, *dpflux; + + assert (W->conn != NULL); + assert (W->ctrl != NULL); + assert (W->data != NULL); + + WI = W->data->WI; + gpot = W->data->gpot; + Ap = W->data->A; + pmobp = W->data->phasemob; + + np2 = np * np; + + pflux = pimpl->compflux_p; + dpflux = pimpl->compflux_deriv_p; + + for (w = i = 0; w < W->conn->number_of_wells; w++) { + pw = wpress[w]; + + for (; i < W->conn->well_connpos[w + 1]; i++, + gpot += np, Ap += np2, pmobp += np, + pflux += np, dpflux += 2 * np) { + + c = W->conn->well_cells[i]; + dp = pw - cpress[c]; + + compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot, + pimpl->flux_work, + pimpl->flux_work + np); + + /* Component flux = Ap * q*/ + matvec(np, np, Ap, pimpl->flux_work , pflux ); + + /* Derivative = Ap * (dq/dp) */ + matmat(np, 2 , Ap, pimpl->flux_work + np, dpflux); + } + } +} + + + static int count_internal_conn(grid_t *G, int c) { @@ -814,6 +866,7 @@ cfs_tpfa_res_assemble(grid_t *G, const double *trans, const double *gravcap_f, const double *cpress, + const double *wpress, const double *porevol, struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ @@ -828,6 +881,11 @@ cfs_tpfa_res_assemble(grid_t *G, compute_compflux_and_deriv(G, cq->nphases, cpress, trans, cq->phasemobf, gravcap_f, cq->Af, h->pimpl); + if ((forces != NULL) && (forces->W != NULL)) { + compute_well_compflux_and_deriv(forces->W, cq->nphases, + cpress, wpress, h->pimpl); + } + np2 = cq->nphases * cq->nphases; for (c = 0; c < G->number_of_cells; c++, zc += cq->nphases) { diff --git a/src/cfs_tpfa_residual.h b/src/cfs_tpfa_residual.h index ac75e455..8e827746 100644 --- a/src/cfs_tpfa_residual.h +++ b/src/cfs_tpfa_residual.h @@ -72,6 +72,7 @@ cfs_tpfa_res_assemble(grid_t *G, const double *trans, const double *gravcap_f, const double *cpress, + const double *wpress, const double *porevol, struct cfs_tpfa_res_data *h); From 3085a4d12b8b72fed8fc68d608bbf8ded464f37d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 23 Nov 2011 23:06:16 +0100 Subject: [PATCH 36/36] Implement csrmatrix_zero() in terms of vector_zero(). --- src/sparse_sys.c | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/sparse_sys.c b/src/sparse_sys.c index 793941bf..7d0583d2 100644 --- a/src/sparse_sys.c +++ b/src/sparse_sys.c @@ -182,9 +182,7 @@ void csrmatrix_zero(struct CSRMatrix *A) /* ---------------------------------------------------------------------- */ { - size_t i; - - for (i = 0; i < A->nnz; i++) { A->sa[i] = 0.0; } + vector_zero(A->nnz, A->sa); }