From c026c35bf75458b372218396d7b9da38a999342c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 13 Oct 2010 18:35:15 +0200 Subject: [PATCH 01/17] Use canonical include guards. Suggested by atgeirr. Template: OPM__HEADER_INCLUDED --- blas_lapack.h | 6 +++--- coarse_conn.h | 6 +++--- coarse_sys.h | 6 +++--- dfs.h | 4 ++-- flow_bc.h | 6 +++--- grid.h | 6 +++--- hash_set.h | 6 +++--- hybsys.h | 6 +++--- hybsys_global.h | 6 +++--- ifsh.h | 6 +++--- ifsh_ms.h | 6 +++--- mimetic.h | 6 +++--- partition.h | 4 ++-- sparse_sys.h | 6 +++--- well.h | 6 +++--- 15 files changed, 43 insertions(+), 43 deletions(-) diff --git a/blas_lapack.h b/blas_lapack.h index d12621ec..ac66f625 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef BLAS_LAPACK_H_INCLUDED -#define BLAS_LAPACK_H_INCLUDED +#ifndef OPM_BLAS_LAPACK_HEADER_INCLUDED +#define OPM_BLAS_LAPACK_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -101,4 +101,4 @@ double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx, } #endif -#endif /* BLAS_LAPACK_H_INCLUDED */ +#endif /* OPM_BLAS_LAPACK_HEADER_INCLUDED */ diff --git a/coarse_conn.h b/coarse_conn.h index 67b375ba..f6ce6349 100644 --- a/coarse_conn.h +++ b/coarse_conn.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef COARSE_CONN_H_INCLUDED -#define COARSE_CONN_H_INCLUDED +#ifndef OPM_COARSE_CONN_HEADER_INCLUDED +#define OPM_COARSE_CONN_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -51,4 +51,4 @@ coarse_topology_destroy(struct coarse_topology *t); } #endif -#endif /* COARSE_CONN_H_INCLUDED */ +#endif /* OPM_COARSE_CONN_HEADER_INCLUDED */ diff --git a/coarse_sys.h b/coarse_sys.h index fc736627..eee06d1d 100644 --- a/coarse_sys.h +++ b/coarse_sys.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef COARSE_SYS_H_INCLUDED -#define COARSE_SYS_H_INCLUDED +#ifndef OPM_COARSE_SYS_HEADER_INCLUDED +#define OPM_COARSE_SYS_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -85,4 +85,4 @@ coarse_sys_compute_Binv(int nb, } #endif -#endif /* COARSE_SYS_H_INCLUDED */ +#endif /* OPM_COARSE_SYS_HEADER_INCLUDED */ diff --git a/dfs.h b/dfs.h index 43747ca5..5213309d 100644 --- a/dfs.h +++ b/dfs.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef DFS_H_INCLUDED -#define DFS_H_INCLUDED +#ifndef OPM_DFS_HEADER_INCLUDED +#define OPM_DFS_HEADER_INCLUDED #ifdef __cplusplus extern "C" { diff --git a/flow_bc.h b/flow_bc.h index 30430a03..6da3e81c 100644 --- a/flow_bc.h +++ b/flow_bc.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef FLOW_BC_H_INCLUDED -#define FLOW_BC_H_INCLUDED +#ifndef OPM_FLOW_BC_HEADER_INCLUDED +#define OPM_FLOW_BC_HEADER_INCLUDED #include @@ -46,4 +46,4 @@ deallocate_flowbc(flowbc_t *fbc); } #endif -#endif /* FLOW_BC_H_INCLUDED */ +#endif /* OPM_FLOW_BC_HEADER_INCLUDED */ diff --git a/grid.h b/grid.h index fd377629..4b575115 100644 --- a/grid.h +++ b/grid.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_GRID_H_INCLUDED -#define OPM_GRID_H_INCLUDED +#ifndef OPM_GRID_HEADER_INCLUDED +#define OPM_GRID_HEADER_INCLUDED #ifdef __cplusplus @@ -66,4 +66,4 @@ typedef struct { } #endif -#endif /* OPM_GRID_H_INCLUDED */ +#endif /* OPM_GRID_HEADER_INCLUDED */ diff --git a/hash_set.h b/hash_set.h index 5a9d8a78..1edfd1ce 100644 --- a/hash_set.h +++ b/hash_set.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef HASH_SET_H_INCLUDED -#define HASH_SET_H_INCLUDED +#ifndef OPM_HASH_SET_HEADER_INCLUDED +#define OPM_HASH_SET_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -67,4 +67,4 @@ hash_set_count_elms(const struct hash_set *set); } #endif -#endif /* HASH_SET_H_INCLUDED */ +#endif /* OPM_HASH_SET_HEADER_INCLUDED */ diff --git a/hybsys.h b/hybsys.h index a476ae03..7e0da807 100644 --- a/hybsys.h +++ b/hybsys.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef HYBSYS_H_INCLUDED -#define HYBSYS_H_INCLUDED +#ifndef OPM_HYBSYS_HEADER_INCLUDED +#define OPM_HYBSYS_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -154,4 +154,4 @@ hybsys_compute_press_flux_well(int nc, const int *pgconn, int nf, } #endif -#endif /* HYBSYS_H_INCLUDED */ +#endif /* OPM_HYBSYS_HEADER_INCLUDED */ diff --git a/hybsys_global.h b/hybsys_global.h index 25ade66b..5680c577 100644 --- a/hybsys_global.h +++ b/hybsys_global.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef HYBSYS_GLOBAL_H_INCLUDED -#define HYBSYS_GLOBAL_H_INCLUDED +#ifndef OPM_HYBSYS_GLOBAL_HEADER_INCLUDED +#define OPM_HYBSYS_GLOBAL_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -55,4 +55,4 @@ hybsys_global_assemble_well_sym(int ngconn_tot, } #endif -#endif /* HYBSYS_GLOBAL_H_INCLUDED */ +#endif /* OPM_HYBSYS_GLOBAL_HEADER_INCLUDED */ diff --git a/ifsh.h b/ifsh.h index 2f0a0dfe..cbefe5fd 100644 --- a/ifsh.h +++ b/ifsh.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef IFSH_H_INCLUDED -#define IFHS_H_INCLUDED +#ifndef OPM_IFSH_HEADER_INCLUDED +#define OPM_IFHS_HEADER_INCLUDED #include "grid.h" #include "well.h" @@ -141,4 +141,4 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h, #endif -#endif /* IFSH_H_INCLUDED */ +#endif /* OPM_IFSH_HEADER_INCLUDED */ diff --git a/ifsh_ms.h b/ifsh_ms.h index 8c9ad9fc..acfe64f2 100644 --- a/ifsh_ms.h +++ b/ifsh_ms.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef IFSH_MS_H_INCLUDED -#define IFSH_MS_H_INCLUDED +#ifndef OPM_IFSH_MS_HEADER_INCLUDED +#define OPM_IFSH_MS_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -68,4 +68,4 @@ ifsh_ms_press_flux(grid_t *G, struct ifsh_ms_data *h, } #endif -#endif /* IFSH_MS_H_INCLUDED */ +#endif /* OPM_IFSH_MS_HEADER_INCLUDED */ diff --git a/mimetic.h b/mimetic.h index 1b73aa1c..659cbb28 100644 --- a/mimetic.h +++ b/mimetic.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef MIMETIC_H_INCLUDED -#define MIMETIC_H_INCLUDED +#ifndef OPM_MIMETIC_HEADER_INCLUDED +#define OPM_MIMETIC_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -85,4 +85,4 @@ void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, } #endif -#endif /* MIMETIC_H_INCLUDED */ +#endif /* OPM_MIMETIC_HEADER_INCLUDED */ diff --git a/partition.h b/partition.h index 7a811529..61336e66 100644 --- a/partition.h +++ b/partition.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef PARTITION_H_INCLUDED -#define PARTITION_H_INCLUDED +#ifndef OPM_PARTITION_HEADER_INCLUDED +#define OPM_PARTITION_HEADER_INCLUDED #ifdef __cplusplus extern "C" { diff --git a/sparse_sys.h b/sparse_sys.h index e9b9d75d..b53491e0 100644 --- a/sparse_sys.h +++ b/sparse_sys.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef SPARSE_SYS_H_INCLUDED -#define SPARSE_SYS_H_INCLUDED +#ifndef OPM_SPARSE_SYS_HEADER_INCLUDED +#define OPM_SPARSE_SYS_HEADER_INCLUDED #include @@ -76,4 +76,4 @@ csrmatrix_zero(struct CSRMatrix *A); } #endif -#endif /* SPARSE_SYS_H_INCLUDED */ +#endif /* OPM_SPARSE_SYS_HEADER_INCLUDED */ diff --git a/well.h b/well.h index f2e3ceea..d4d0dae9 100644 --- a/well.h +++ b/well.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef WELL_H_INCLUDED -#define WELL_H_INCLUDED +#ifndef OPM_WELL_HEADER_INCLUDED +#define OPM_WELL_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -54,4 +54,4 @@ derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells); } #endif -#endif /* WELL_H_INCLUDED */ +#endif /* OPM_WELL_HEADER_INCLUDED */ From fdf089c6a8d9e8de1e0c135ccdf9d516415f75ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 13 Oct 2010 19:02:06 +0200 Subject: [PATCH 02/17] Remark that enumerate_active_bf() is restricted to internal CFs. --- coarse_sys.c | 1 + 1 file changed, 1 insertion(+) diff --git a/coarse_sys.c b/coarse_sys.c index 793e632d..d30f4fbc 100644 --- a/coarse_sys.c +++ b/coarse_sys.c @@ -311,6 +311,7 @@ enumerate_active_bf(struct coarse_topology *ct, b_out = (b1 == b_in) ? b2 : b1; if (b_out >= 0) { + /* Restricted to internal cf's for now. */ m->bfno[cf] = act++; } } From 2985200ad58362800c4ea7e0b0a54a94704002d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 13 Oct 2010 23:21:17 +0200 Subject: [PATCH 03/17] Ignore files generated by the Autotools. --- .hgignore | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 .hgignore diff --git a/.hgignore b/.hgignore new file mode 100644 index 00000000..b9c7ed97 --- /dev/null +++ b/.hgignore @@ -0,0 +1,12 @@ +syntax: glob +**Makefile.in +aclocal.m4 +autom4te.cache +config.* +configure +depcomp +install-sh +ltmain.sh +m4/libtool.m4 +m4/lt*.m4 +missing From 4a2b5e4e0eb9a34ab8fe6623520122837628ab49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 14 Oct 2010 13:58:47 +0200 Subject: [PATCH 04/17] Move implementation of vector_zero() to central location. There is no need for each file to contain a separate, though trivial, implementation of this feature. --- coarse_sys.c | 13 ------------- ifsh_ms.c | 11 ----------- sparse_sys.c | 13 +++++++++++++ sparse_sys.h | 6 ++++++ 4 files changed, 19 insertions(+), 24 deletions(-) diff --git a/coarse_sys.c b/coarse_sys.c index d30f4fbc..780cffd8 100644 --- a/coarse_sys.c +++ b/coarse_sys.c @@ -1054,19 +1054,6 @@ define_csr_sparsity(size_t nc, size_t m, struct bf_asm_data *bf_asm) } -/* ---------------------------------------------------------------------- */ -/* v = zeros([n, 1]) */ -/* ---------------------------------------------------------------------- */ -static void -vector_zero(size_t n, double *v) -/* ---------------------------------------------------------------------- */ -{ - size_t i; - - for (i = 0; i < n; i++) { v[i] = 0.0; } -} - - /* ---------------------------------------------------------------------- */ /* Assemble system of linear equations corresponding to local * discretisation of flow problem on domain connected to coarse face diff --git a/ifsh_ms.c b/ifsh_ms.c index 0a6a5217..a4ac7335 100644 --- a/ifsh_ms.c +++ b/ifsh_ms.c @@ -348,17 +348,6 @@ ifsh_ms_impl_construct(grid_t *G , } -/* ---------------------------------------------------------------------- */ -static void -vector_zero(size_t m, double *v) -/* ---------------------------------------------------------------------- */ -{ - size_t i; - - for (i = 0; i < m; i++) { v[i] = 0.0; } -} - - /* ---------------------------------------------------------------------- */ static void average_flux(size_t nf, const int *N, double *flux) /* Poor name */ diff --git a/sparse_sys.c b/sparse_sys.c index 7d48f06d..8fbfed4b 100644 --- a/sparse_sys.c +++ b/sparse_sys.c @@ -192,3 +192,16 @@ csrmatrix_zero(struct CSRMatrix *A) for (i = 0; i < A->nnz; i++) { A->sa[i] = 0.0; } } + + +/* ---------------------------------------------------------------------- */ +/* v = zeros([n, 1]) */ +/* ---------------------------------------------------------------------- */ +void +vector_zero(size_t n, double *v) +/* ---------------------------------------------------------------------- */ +{ + size_t i; + + for (i = 0; i < n; i++) { v[i] = 0.0; } +} diff --git a/sparse_sys.h b/sparse_sys.h index b53491e0..aadaedf5 100644 --- a/sparse_sys.h +++ b/sparse_sys.h @@ -72,6 +72,12 @@ csrmatrix_delete(struct CSRMatrix *A); void csrmatrix_zero(struct CSRMatrix *A); +/* ---------------------------------------------------------------------- */ +/* v = zeros([n, 1]) */ +/* ---------------------------------------------------------------------- */ +void +vector_zero(size_t n, double *v); + #ifdef __cplusplus } #endif From 84a6a1c9aef28da73ee4a51d6985033f117d75aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 14 Oct 2010 14:20:11 +0200 Subject: [PATCH 05/17] Add facilities for computing/updating gpress/Binv. We were already computing the inverse IP, but now centralise the gpress as well. Moreover, the mim_ip_*_update() functions will assist in the pressure solvers accepting effective inner products and gravity pressures only. --- mimetic.c | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ mimetic.h | 16 ++++++++++++++ 2 files changed, 81 insertions(+) diff --git a/mimetic.c b/mimetic.c index 48e3a0d3..f9b311b8 100644 --- a/mimetic.c +++ b/mimetic.c @@ -183,3 +183,68 @@ mim_ip_linpress_exact(int nf, int nconn, int d, dgemm_("No Transpose", "Transpose", &m, &m, &n, &a1, work, &ld1, N, &ld1, &a2, Binv, &ldBinv); } + + +/* ---------------------------------------------------------------------- */ +void +mim_ip_compute_gpress(int nc, int d, const double *grav, + const int *pconn, const int *conn, + const double *fcentroid, const double *ccentroid, + double *gpress) +/* ---------------------------------------------------------------------- */ +{ + int c, i, j; + + const double *cc, *fc; + + for (c = i = 0; c < nc; c++) { + cc = ccentroid + (c * d); + + for (; i < pconn[c + 1]; i++) { + fc = fcentroid + (conn[i] * d); + + gpress[i] = 0.0; + for (j = 0; j < d; j++) { + gpress[i] += grav[j] * (fc[j] - cc[j]); + } + } + } +} + + +/* inv(B) <- \lambda_t(s)*inv(B) */ +/* ---------------------------------------------------------------------- */ +void +mim_ip_mobility_update(int nc, const int *pconn, const double *totmob, + double *Binv) +/* ---------------------------------------------------------------------- */ +{ + int c, i, n, p2; + + for (c = p2 = 0; c < nc; c++) { + n = pconn[c + 1] - pconn[c]; + + for (i = 0; i < n * n; i++) { + Binv[p2 + i] *= totmob[c]; + } + + p2 += n * n; + } +} + + +/* G <- \sum_i \rho_i f_i(s) * G */ +/* ---------------------------------------------------------------------- */ +void +mim_ip_density_update(int nc, const int *pconn, const double *omega, + double *gpress) +/* ---------------------------------------------------------------------- */ +{ + int c, i; + + for (c = i = 0; c < nc; c++) { + for (; i < pconn[c + 1]; i++) { + gpress[i] *= omega[c]; + } + } +} diff --git a/mimetic.h b/mimetic.h index 659cbb28..82e347cf 100644 --- a/mimetic.h +++ b/mimetic.h @@ -81,6 +81,22 @@ void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, double *farea, double *ccentroid, double *cvol, double *perm, double *Binv); +void +mim_ip_compute_gpress(int nc, int d, const double *grav, + const int *pconn, const int *conn, + const double *fcentroid, const double *ccentroid, + double *gpress); + +/* inv(B) <- \lambda_t(s)*inv(B) */ +void +mim_ip_mobility_update(int nc, const int *pconn, const double *totmob, + double *Binv); + +/* G <- \sum_i \rho_i f_i(s) * G */ +void +mim_ip_density_update(int nc, const int *pconn, const double *omega, + double *gpress); + #ifdef __cplusplus } #endif From d2fa5f4e103d6e78889493831640031d5b812460 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 15 Oct 2010 18:49:46 +0200 Subject: [PATCH 06/17] Re-factor IFSH into FSH common and IFSH special parts. In preparation of adding compressible flow solver. --- Makefile.am | 2 + fsh_common.c | 237 +++++++++++++++++++++ fsh_common.h | 61 ++++++ fsh_common_impl.h | 79 +++++++ ifsh.c | 522 ++++++++++++++-------------------------------- ifsh.h | 63 ++---- 6 files changed, 556 insertions(+), 408 deletions(-) create mode 100644 fsh_common.c create mode 100644 fsh_common.h create mode 100644 fsh_common_impl.h diff --git a/Makefile.am b/Makefile.am index 40748115..9fb8f72c 100644 --- a/Makefile.am +++ b/Makefile.am @@ -14,6 +14,7 @@ coarse_conn.h \ coarse_sys.h \ dfs.h \ flow_bc.h \ +fsh_common.h \ grid.h \ hash_set.h \ hybsys_global.h \ @@ -33,6 +34,7 @@ coarse_conn.c \ coarse_sys.c \ dfs.c \ flow_bc.c \ +fsh_common.c \ hash_set.c \ hybsys.c \ hybsys_global.c \ diff --git a/fsh_common.c b/fsh_common.c new file mode 100644 index 00000000..8561b17d --- /dev/null +++ b/fsh_common.c @@ -0,0 +1,237 @@ +/* + 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 +#include +#include +#include +#include + +#include "grid.h" +#include "flow_bc.h" +#include "well.h" + +#include "fsh_common.h" +#include "fsh_common_impl.h" + +#include "hybsys.h" +#include "hybsys_global.h" + +#if defined MAX +#undef MAX +#endif +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + + + +/* ---------------------------------------------------------------------- */ +/* Release dynamic memory resources associated to internal data of a + * particular (incompressible) flow solver instance. */ +/* ---------------------------------------------------------------------- */ +static void +fsh_destroy_impl(struct fsh_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + if (pimpl != NULL) { + hybsys_well_free(pimpl->wsys ); + hybsys_free (pimpl->sys ); + free (pimpl->ddata); + free (pimpl->idata); + } + + free(pimpl); +} + + +/* ---------------------------------------------------------------------- */ +/* Eliminate 'npp' prescribed (Dirichlet) conditions (locdof,dofval) + * from global system of linear equations. Move known values to RHS + * whilst zeroing coefficient matrix contributions. Local system of + * dimension 'n'. */ +/* ---------------------------------------------------------------------- */ +static void +fsh_explicit_elimination(int n, int npp, + int *locdof, double *dofval, + double *S, double *r) +/* ---------------------------------------------------------------------- */ +{ + int i, p; + + for (i = 0; i < npp; i++) { + for (p = (locdof[i] + 1) % n; p != locdof[i]; p = (p + 1) % n) { + /* Subtract known values from RHS. */ + r[p] -= S[p + locdof[i]*n] * dofval[i]; + + /* Eliminate DOF (zero row/col "locdof[i]"; leave diagonal). */ + S[p + locdof[i]*n] = 0.0; + S[locdof[i] + p*n] = 0.0; + } + + /* Produce trivial equation S(i,i)*x(i) = S(i,i)*x0(i). */ + r[p] = S[p * (n + 1)] * dofval[i]; + } +} + + +/* ---------------------------------------------------------------------- */ +/* Release memory resources for ifsh data-handle 'h' */ +/* ---------------------------------------------------------------------- */ +void +fsh_destroy(struct fsh_data *h) +/* ---------------------------------------------------------------------- */ +{ + if (h != NULL) { + fsh_destroy_impl(h->pimpl); + csrmatrix_delete(h->A ); + } + + free(h); +} + + +/* ---------------------------------------------------------------------- */ +struct fsh_impl * +fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz) +/* ---------------------------------------------------------------------- */ +{ + struct fsh_impl *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->idata = malloc(idata_sz * sizeof *new->idata); + new->ddata = malloc(ddata_sz * sizeof *new->ddata); + + new->sys = NULL; + new->wsys = NULL; + + if ((new->idata == NULL) || (new->ddata == NULL)) { + fsh_destroy_impl(new); + new = NULL; + } + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +/* Determine nnz (=sum(diff(facePos)^2)) and max(diff(facePos) for grid */ +/* ---------------------------------------------------------------------- */ +void +fsh_count_grid_dof(grid_t *G, int *max_ngdof, size_t *sum_ngdof2) +/* ---------------------------------------------------------------------- */ +{ + int c, n; + + *max_ngdof = INT_MIN; + *sum_ngdof2 = 0; + + for (c = 0; c < G->number_of_cells; c++) { + n = G->cell_facepos[c + 1] - G->cell_facepos[c]; + + *max_ngdof = MAX(*max_ngdof, n); + *sum_ngdof2 += n * n; + } +} + + +/* ---------------------------------------------------------------------- */ +/* Impose boundary conditions on local contribution to global system. */ +/* ---------------------------------------------------------------------- */ +int +fsh_impose_bc(int nconn, int *conn, flowbc_t *bc, + struct fsh_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + int i, npp, f; + + npp = 0; + for (i = 0; i < nconn; i++) { + f = conn[i]; + + if (bc->type[f] == PRESSURE) { + pimpl->work [npp] = bc->bcval[f]; + pimpl->iwork[npp] = i; + + npp += 1; + } else if (bc->type[f] == FLUX) { + pimpl->sys->r[i] -= bc->bcval[f]; + } + } + + if (npp > 0) { + fsh_explicit_elimination(nconn, npp, + pimpl->iwork, + pimpl->work, + pimpl->sys->S, + pimpl->sys->r); + } + + return npp; +} + + +/* ---------------------------------------------------------------------- */ +void +fsh_define_impl_arrays(size_t nc, + size_t nnu, + size_t nhf, + size_t max_ncf, + well_t *W, + struct fsh_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + pimpl->cflux = pimpl->ddata + 2 * nnu; + pimpl->work = pimpl->cflux + nhf; + + pimpl->gdof_pos = pimpl->idata; + pimpl->gdof = pimpl->gdof_pos + (nc + 1); + pimpl->iwork = pimpl->gdof + nhf; + + if (W != NULL) { + pimpl->cwell_pos = pimpl->iwork + max_ncf; + pimpl->cwells = pimpl->cwell_pos + nc + 1; + + pimpl->WI = pimpl->work + max_ncf; + pimpl->wdp = pimpl->WI + W->well_connpos[ W->number_of_wells ]; + } +} + + +/* ---------------------------------------------------------------------- */ +void +fsh_define_cell_wells(size_t nc, well_t *W, struct fsh_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + memset(pimpl->cwell_pos, 0, (nc + 1) * sizeof *pimpl->cwell_pos); + + derive_cell_wells(nc, W, pimpl->cwell_pos, pimpl->cwells); +} + + +/* ---------------------------------------------------------------------- */ +void +fsh_define_linsys_arrays(struct fsh_data *h) +/* ---------------------------------------------------------------------- */ +{ + h->b = h->pimpl->ddata; + h->x = h->b + h->A->m; +} diff --git a/fsh_common.h b/fsh_common.h new file mode 100644 index 00000000..9ba62402 --- /dev/null +++ b/fsh_common.h @@ -0,0 +1,61 @@ +/* + 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_IFSH_HEADER_INCLUDED +#define OPM_IFHS_HEADER_INCLUDED + +#ifdef __cplusplus +extern "C" { +#endif + +struct CSRMatrix; +struct fsh_impl; + +/** Contains the linear system for assembly, as well as internal data + * for the assembly routines. + */ +struct fsh_data { + /* Let \f$n_i\f$ be the number of connections/faces of grid cell + * number \f$i\f$. Then max_ngconn = \f$\max_i n_i\f$ + */ + int max_ngconn; + /* With n_i as above, sum_ngconn2 = \f$\sum_i n_i^2\f$ */ + size_t sum_ngconn2; + + /* Linear system */ + struct CSRMatrix *A; /* Coefficient matrix */ + double *b; /* System RHS */ + double *x; /* Solution */ + + /* Private implementational details. */ + struct fsh_impl *pimpl; +}; + + + +/** Destroys the fsh data object */ +void +fsh_destroy(struct fsh_data *h); + +#ifdef __cplusplus +} +#endif + + +#endif /* OPM_IFSH_HEADER_INCLUDED */ diff --git a/fsh_common_impl.h b/fsh_common_impl.h new file mode 100644 index 00000000..04992770 --- /dev/null +++ b/fsh_common_impl.h @@ -0,0 +1,79 @@ +/* + 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_FSH_COMMON_IMPL_HEADER_INCLUDED +#define OPM_FSH_COMMON_IMPL_HEADER_INCLUDED + +/* Internal header. Don't install. */ + +struct fsh_impl { + int nc, nf, nw; /* Number of cells, faces, wells */ + + /* Topology */ + int *gdof_pos; /* Pointers, grid DOFs (== cell_facepos) */ + int *gdof; /* Grid DOFs (== cell_faces) */ + + int *cwell_pos; /* Start pointers, well DOFs (c->w) */ + int *cwells; /* Well DOFs (c->w) */ + + /* Discretisation data */ + double *WI; /* Permuted well production indices */ + double *wdp; /* Permuted well gravity pressures */ + + double *cflux; /* Cell (half-contact) fluxes */ + + struct hybsys *sys; /* Hybrid cell contribs */ + struct hybsys_well *wsys; /* Hybrid cell contribs from wells */ + + double *work; /* Scratch array, floating point */ + int *iwork; /* Scratch array, integers */ + + /* Linear storage goes here... */ + int *idata; /* Actual storage array, integers */ + double *ddata; /* Actual storage array, floating point */ +}; + + +struct fsh_impl * +fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz); + +void +fsh_count_grid_dof(grid_t *G, int *max_ngdof, size_t *sum_ngdof2); + +int +fsh_impose_bc(int ndof, + int *dof, + flowbc_t *bc, + struct fsh_impl *pimpl); + +void +fsh_define_impl_arrays(size_t nc, + size_t nnu, + size_t nhf, + size_t max_ncf, + well_t *W, + struct fsh_impl *pimpl); + +void +fsh_define_cell_wells(size_t nc, well_t *W, struct fsh_impl *pimpl); + +void +fsh_define_sys_arrays(struct fsh_data *h); + +#endif /* OPM_FSH_COMMON_IMPL_HEADER_INCLUDED */ diff --git a/ifsh.c b/ifsh.c index 900d33c6..2e9bcdea 100644 --- a/ifsh.c +++ b/ifsh.c @@ -24,7 +24,9 @@ #include #include +#include "fsh_common.h" #include "ifsh.h" +#include "fsh_common_impl.h" #include "hybsys.h" #include "hybsys_global.h" @@ -33,59 +35,10 @@ #endif #define MAX(a,b) (((a) > (b)) ? (a) : (b)) -struct ifsh_impl { - int nc, nf, nw; /* Number of cells, faces, wells */ - - double *Binv; /* Effective inverse IP */ - double *gpress; /* Effective gravity pressure */ - double *cflux; /* Cell (half-contact) fluxes */ - - double *WI; /* Effective inverse IP, wells */ - double *wdp; /* Effective grav. press, wells */ - - int *pgconn; /* Start pointers, grid connections */ - int *gconn; /* Grid connections */ - - int *cwpos; /* Start pointers, c->w mapping */ - int *cwells; /* c->w mapping */ - - struct hybsys *sys; /* Hybrid cell contribs */ - struct hybsys_well *wsys; /* Hybrid cell contribs from wells */ - - double *work; /* Scratch array, floating point */ - int *iwork; /* Scratch array, integers */ - - /* Linear storage goes here... */ - int *idata; /* Actual storage array, integers */ - double *ddata; /* Actual storage array, floating point */ -}; - - -/* ---------------------------------------------------------------------- */ -/* Determine nnz (=sum(diff(facePos)^2)) and max(diff(facePos) for grid */ -/* ---------------------------------------------------------------------- */ -static void -count_grid_connections(grid_t *G, int *max_ngconn, size_t *sum_ngconn2) -/* ---------------------------------------------------------------------- */ -{ - int c, n; - - *max_ngconn = INT_MIN; - *sum_ngconn2 = 0; - - for (c = 0; c < G->number_of_cells; c++) { - n = G->cell_facepos[c + 1] - G->cell_facepos[c]; - - *max_ngconn = MAX(*max_ngconn, n); - *sum_ngconn2 += n * n; - } -} - /* ---------------------------------------------------------------------- */ static void -ifsh_compute_table_sz(grid_t *G, well_t *W, - int max_ngconn, size_t sum_ngconn2, +ifsh_compute_table_sz(grid_t *G, well_t *W, int max_ngconn, size_t *nnu, size_t *idata_sz, size_t *ddata_sz) /* ---------------------------------------------------------------------- */ { @@ -96,18 +49,18 @@ ifsh_compute_table_sz(grid_t *G, well_t *W, nc = G->number_of_cells; ngconn_tot = G->cell_facepos[nc]; - *idata_sz = max_ngconn; /* iwork */ + *idata_sz = nc + 1; /* gdof_pos */ + *idata_sz += ngconn_tot; /* gdof */ + *idata_sz += max_ngconn; /* iwork */ *ddata_sz = 2 * (*nnu); /* rhs + soln */ - *ddata_sz += sum_ngconn2; /* Binv */ - *ddata_sz += ngconn_tot; /* gpress */ *ddata_sz += ngconn_tot; /* cflux */ *ddata_sz += max_ngconn; /* work */ if (W != NULL) { *nnu += W->number_of_wells; - /* cwpos */ + /* cwell_pos */ *idata_sz += nc + 1; /* cwells */ @@ -122,277 +75,59 @@ ifsh_compute_table_sz(grid_t *G, well_t *W, } -/* ---------------------------------------------------------------------- */ -/* Allocate and define supporting structures for assembling the global - * system of linear equations to couple the grid (reservoir) - * connections represented by 'G' and, if present (i.e., non-NULL), - * the well connections represented by 'W'. */ -/* ---------------------------------------------------------------------- */ -struct ifsh_data * -ifsh_construct(grid_t *G, well_t *W) -/* ---------------------------------------------------------------------- */ -{ - int nc, ngconn_tot; - size_t idata_sz, ddata_sz, nnu; - struct ifsh_data *new; - - assert (G != NULL); - - new = malloc(1 * sizeof *new); - if (new != NULL) { - new->A = hybsys_define_globconn(G, W); - new->pimpl = NULL; - - if (new->A == NULL) { - ifsh_destroy(new); - new = NULL; - } - } - - if (new != NULL) { - new->pimpl = malloc(1 * sizeof *new->pimpl); - - if (new->pimpl == NULL) { - ifsh_destroy(new); - new = NULL; - } else { - new->pimpl->sys = NULL; new->pimpl->wsys = NULL; - new->pimpl->idata = NULL; new->pimpl->ddata = NULL; - } - } - - if (new != NULL) { - count_grid_connections(G, &new->max_ngconn, &new->sum_ngconn2); - - ifsh_compute_table_sz(G, W, new->max_ngconn, new->sum_ngconn2, - &nnu, &idata_sz, &ddata_sz); - - nc = G->number_of_cells; - ngconn_tot = G->cell_facepos[nc]; - - new->pimpl->idata = malloc(idata_sz * sizeof *new->pimpl->idata); - new->pimpl->ddata = malloc(ddata_sz * sizeof *new->pimpl->ddata); - new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn, - nc, ngconn_tot); - if ((new->pimpl->idata == NULL) || - (new->pimpl->ddata == NULL) || (new->pimpl->sys == NULL)) { - ifsh_destroy(new); - new = NULL; - } else { - new->b = new->pimpl->ddata; - new->x = new->b + nnu; - new->pimpl->Binv = new->x + nnu; - new->pimpl->gpress = new->pimpl->Binv + new->sum_ngconn2; - new->pimpl->cflux = new->pimpl->gpress + ngconn_tot; - new->pimpl->work = new->pimpl->cflux + ngconn_tot; - - new->pimpl->iwork = new->pimpl->idata; - - hybsys_init(new->max_ngconn, new->pimpl->sys); - } - } - - if ((new != NULL) && (W != NULL)) { - new->pimpl->cwpos = new->pimpl->iwork + new->max_ngconn; - new->pimpl->cwells = new->pimpl->cwpos + nc + 1; - - memset(new->pimpl->cwpos, 0, (nc + 1) * sizeof *new->pimpl->cwpos); - - derive_cell_wells(nc, W, new->pimpl->cwpos, new->pimpl->cwells); - - new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc, - new->pimpl->cwpos); - - if (new->pimpl->wsys == NULL) { - ifsh_destroy(new); - new = NULL; - } else { - new->pimpl->WI = new->pimpl->work + new->max_ngconn; - new->pimpl->wdp = new->pimpl->WI + W->well_connpos[ W->number_of_wells ]; - } - } - - if (new != NULL) { - new->pimpl->nc = nc; - new->pimpl->nf = G->number_of_faces; - new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0; - - /* Contentious. Imposes severe restrictions on lifetime(G). */ - new->pimpl->pgconn = G->cell_facepos; - new->pimpl->gconn = G->cell_faces; - } - - return new; -} - - -/* ---------------------------------------------------------------------- */ -/* Release dynamic memory resources associated to internal data of a - * particular (incompressible) flow solver instance. */ /* ---------------------------------------------------------------------- */ static void -ifsh_destroy_impl(struct ifsh_impl *pimpl) +ifsh_set_effective_well_params(const double *WI, + const double *wdp, + struct fsh_data *ifsh) /* ---------------------------------------------------------------------- */ { - if (pimpl != NULL) { - hybsys_well_free(pimpl->wsys ); - hybsys_free (pimpl->sys ); - free (pimpl->ddata); - free (pimpl->idata); - } + int c, nc, i, perf; + int *cwpos, *cwells; + double *wsys_WI, *wsys_wdp; - free(pimpl); -} + nc = ifsh->pimpl->nc; + cwpos = ifsh->pimpl->cwell_pos; + cwells = ifsh->pimpl->cwells; + wsys_WI = ifsh->pimpl->WI; + wsys_wdp = ifsh->pimpl->wdp; + for (c = i = 0; c < nc; c++) { + for (; i < cwpos[c + 1]; i++) { + perf = cwells[2*i + 1]; -/* ---------------------------------------------------------------------- */ -/* Release memory resources for ifsh data-handle 'h' */ -/* ---------------------------------------------------------------------- */ -void -ifsh_destroy(struct ifsh_data *h) -/* ---------------------------------------------------------------------- */ -{ - if (h != NULL) { - ifsh_destroy_impl(h->pimpl); - csrmatrix_delete (h->A ); - } - - free(h); -} - - -/* ---------------------------------------------------------------------- */ -/* Eliminate 'npp' prescribed (Dirichlet) conditions (locdof,dofval) - * from global system of linear equations. Move known values to RHS - * whilst zeroing coefficient matrix contributions. Local system of - * dimension 'n'. */ -/* ---------------------------------------------------------------------- */ -static void -ifsh_explicit_elimination(int n, int npp, - int *locdof, double *dofval, - double *S, double *r) -/* ---------------------------------------------------------------------- */ -{ - int i, p; - - for (i = 0; i < npp; i++) { - for (p = (locdof[i] + 1) % n; p != locdof[i]; p = (p + 1) % n) { - /* Subtract known values from RHS. */ - r[p] -= S[p + locdof[i]*n] * dofval[i]; - - /* Eliminate DOF (zero row/col "locdof[i]"; leave diagonal). */ - S[p + locdof[i]*n] = 0.0; - S[locdof[i] + p*n] = 0.0; - } - - /* Produce trivial equation S(i,i)*x(i) = S(i,i)*x0(i). */ - r[p] = S[p * (n + 1)] * dofval[i]; - } -} - - -/* ---------------------------------------------------------------------- */ -/* Impose boundary conditions on local contribution to global system. */ -/* ---------------------------------------------------------------------- */ -static int -ifsh_impose_bc(int nconn, int *conn, flowbc_t *bc, - struct ifsh_impl *pimpl) -/* ---------------------------------------------------------------------- */ -{ - int i, npp, f; - - npp = 0; - for (i = 0; i < nconn; i++) { - f = conn[i]; - - if (bc->type[f] == PRESSURE) { - pimpl->work [npp] = bc->bcval[f]; - pimpl->iwork[npp] = i; - - npp += 1; - } else if (bc->type[f] == FLUX) { - pimpl->sys->r[i] -= bc->bcval[f]; + wsys_WI [i] = WI [perf]; + wsys_wdp[i] = wdp[perf]; } } - - if (npp > 0) { - ifsh_explicit_elimination(nconn, npp, - pimpl->iwork, - pimpl->work, - pimpl->sys->S, - pimpl->sys->r); - } - - return npp; } -/* ---------------------------------------------------------------------- */ -static void -ifsh_set_effective_params(const double *Binv, - const double *gpress, - const double *totmob, - const double *omega, - struct ifsh_data *ifsh) -/* ---------------------------------------------------------------------- */ -{ - int c, n, nc, p1, p2, i; - int *pgconn, *gconn; - double *BI, *gp; - - nc = ifsh->pimpl->nc; - pgconn = ifsh->pimpl->pgconn; - gconn = ifsh->pimpl->gconn; - - BI = ifsh->pimpl->Binv; - gp = ifsh->pimpl->gpress; - - p1 = p2 = 0; - for (c = 0; c < nc; c++) { - n = pgconn[c + 1] - pgconn[c]; - - for (i = 0; i < n ; i++) { - gp[p1 + i] = omega[c] * gpress[p1 + i]; - } - - for (i = 0; i < n * n; i++) { - BI[p2 + i] = totmob[c] * Binv[p2 + i]; - } - - p1 += n; - p2 += n * n; - } -} - - - /* ---------------------------------------------------------------------- */ static int -ifsh_assemble_grid(flowbc_t *bc, - const double *src, - struct ifsh_data *ifsh) +ifsh_assemble_grid(flowbc_t *bc, + const double *Binv, + const double *gpress, + const double *src, + struct fsh_data *ifsh) /* ---------------------------------------------------------------------- */ { int c, n, nc, p1, p2; int npp; int *pgconn, *gconn; - double *BI, *gp; nc = ifsh->pimpl->nc; - pgconn = ifsh->pimpl->pgconn; - gconn = ifsh->pimpl->gconn; - - BI = ifsh->pimpl->Binv; - gp = ifsh->pimpl->gpress; + pgconn = ifsh->pimpl->gdof_pos; + gconn = ifsh->pimpl->gdof; p1 = p2 = npp = 0; for (c = 0; c < nc; c++) { n = pgconn[c + 1] - pgconn[c]; - hybsys_cellcontrib_symm(c, n, p1, p2, gp, src, BI, + hybsys_cellcontrib_symm(c, n, p1, p2, gpress, src, Binv, ifsh->pimpl->sys); - npp += ifsh_impose_bc(n, gconn + p1, bc, ifsh->pimpl); + npp += fsh_impose_bc(n, gconn + p1, bc, ifsh->pimpl); hybsys_global_assemble_cell(n, gconn + p1, ifsh->pimpl->sys->S, @@ -409,40 +144,10 @@ ifsh_assemble_grid(flowbc_t *bc, /* ---------------------------------------------------------------------- */ static void -ifsh_set_effective_well_params(const double *WI, - const double *wdp, - const double *totmob, - const double *omega, - struct ifsh_data *ifsh) -/* ---------------------------------------------------------------------- */ -{ - int c, nc, i, perf; - int *cwpos, *cwells; - double *wsys_WI, *wsys_wdp; - - nc = ifsh->pimpl->nc; - cwpos = ifsh->pimpl->cwpos; - cwells = ifsh->pimpl->cwells; - wsys_WI = ifsh->pimpl->WI; - wsys_wdp = ifsh->pimpl->wdp; - - for (c = i = 0; c < nc; c++) { - for (; i < cwpos[c + 1]; i++) { - perf = cwells[2*i + 1]; - - wsys_WI [i] = totmob[c] * WI [perf]; - wsys_wdp[i] = omega [c] * wdp[perf]; - } - } -} - - -/* ---------------------------------------------------------------------- */ -static void -ifsh_impose_well_control(int c, - flowbc_t *bc, - well_control_t *wctrl, - struct ifsh_data *ifsh) +ifsh_impose_well_control(int c, + flowbc_t *bc, + well_control_t *wctrl, + struct fsh_data *ifsh) /* ---------------------------------------------------------------------- */ { int ngconn, nwconn, i, w1, w2, wg, f; @@ -454,10 +159,10 @@ ifsh_impose_well_control(int c, /* Enforce symmetric system */ assert (ifsh->pimpl->wsys->r2w == ifsh->pimpl->wsys->w2r); - pgconn = ifsh->pimpl->pgconn; - pwconn = ifsh->pimpl->cwpos; + pgconn = ifsh->pimpl->gdof_pos; + pwconn = ifsh->pimpl->cwell_pos; - gconn = ifsh->pimpl->gconn + pgconn[c]; + gconn = ifsh->pimpl->gdof + pgconn[c]; wconn = ifsh->pimpl->cwells + 2*pwconn[c]; ngconn = pgconn[c + 1] - pgconn[c]; @@ -517,9 +222,9 @@ ifsh_impose_well_control(int c, /* ---------------------------------------------------------------------- */ static int -ifsh_assemble_well(flowbc_t *bc, - well_control_t *wctrl, - struct ifsh_data *ifsh) +ifsh_assemble_well(flowbc_t *bc, + well_control_t *wctrl, + struct fsh_data *ifsh) /* ---------------------------------------------------------------------- */ { int npp; @@ -529,9 +234,9 @@ ifsh_assemble_well(flowbc_t *bc, nc = ifsh->pimpl->nc; - pgconn = ifsh->pimpl->pgconn; - gconn = ifsh->pimpl->gconn; - pwconn = ifsh->pimpl->cwpos; + pgconn = ifsh->pimpl->gdof_pos; + gconn = ifsh->pimpl->gdof; + pwconn = ifsh->pimpl->cwell_pos; wconn = ifsh->pimpl->cwells; for (c = 0; c < nc; c++) { @@ -575,6 +280,105 @@ ifsh_assemble_well(flowbc_t *bc, } +/* ====================================================================== + * Public routines follow. + * ====================================================================== */ + + +/* ---------------------------------------------------------------------- */ +/* Allocate and define supporting structures for assembling the global + * system of linear equations to couple the grid (reservoir) + * connections represented by 'G' and, if present (i.e., non-NULL), + * the well connections represented by 'W'. */ +/* ---------------------------------------------------------------------- */ +struct fsh_data * +ifsh_construct(grid_t *G, well_t *W) +/* ---------------------------------------------------------------------- */ +{ + int nc, ngconn_tot; + size_t idata_sz, ddata_sz, nnu; + struct fsh_data *new; + + assert (G != NULL); + + /* Allocate master structure, define system matrix sparsity */ + new = malloc(1 * sizeof *new); + if (new != NULL) { + new->A = hybsys_define_globconn(G, W); + new->pimpl = NULL; + + if (new->A == NULL) { + fsh_destroy(new); + new = NULL; + } + } + + + /* Allocate implementation structure */ + if (new != NULL) { + fsh_count_grid_dof(G, &new->max_ngconn, &new->sum_ngconn2); + + ifsh_compute_table_sz(G, W, new->max_ngconn, + &nnu, &idata_sz, &ddata_sz); + + new->pimpl = fsh_impl_allocate_basic(idata_sz, ddata_sz); + + if (new->pimpl == NULL) { + fsh_destroy(new); + new = NULL; + } + } + + + /* Allocate Schur complement contributions. Symmetric system. */ + if (new != NULL) { + nc = G->number_of_cells; + ngconn_tot = G->cell_facepos[nc]; + + fsh_define_linsys_arrays(new); + fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn, + W, new->pimpl); + + new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn, + nc, ngconn_tot); + + if (W != NULL) { + fsh_define_cell_wells(nc, W, new->pimpl); + + new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc, + new->pimpl->cwell_pos); + } + + if ((new->pimpl->sys == NULL) || + ((W != NULL) && (new->pimpl->wsys == NULL))) { + /* Failed to allocate ->sys or ->wsys (if W != NULL) */ + fsh_destroy(new); + new = NULL; + } + } + + + if (new != NULL) { + /* All allocations succeded. Fill metadata and return. */ + new->pimpl->nc = nc; + new->pimpl->nf = G->number_of_faces; + new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0; + + memcpy(new->pimpl->gdof_pos, + G->cell_facepos , + (nc + 1) * sizeof *new->pimpl->gdof_pos); + + memcpy(new->pimpl->gdof , + G->cell_faces , + ngconn_tot * sizeof *new->pimpl->gdof); + + hybsys_init(new->max_ngconn, new->pimpl->sys); + } + + return new; +} + + /* ---------------------------------------------------------------------- */ /* Assemble global system of linear equations * @@ -586,37 +390,32 @@ ifsh_assemble_well(flowbc_t *bc, /* ---------------------------------------------------------------------- */ void ifsh_assemble(flowbc_t *bc, - double *src, - double *Binv, - double *gpress, + const double *src, + const double *Binv, + const double *gpress, well_control_t *wctrl, - double *WI, - double *wdp, - double *totmob, /* \sum_i \lambda_i */ - double *omega, /* \sum_i \rho_i f_i */ - struct ifsh_data *ifsh) + const double *WI, + const double *wdp, + struct fsh_data *ifsh) /* ---------------------------------------------------------------------- */ { int npp; /* Number of prescribed pressure values */ - ifsh_set_effective_params(Binv, gpress, totmob, omega, ifsh); - hybsys_schur_comp_symm(ifsh->pimpl->nc, - ifsh->pimpl->pgconn, - ifsh->pimpl->Binv, - ifsh->pimpl->sys); + ifsh->pimpl->gdof_pos, + Binv, ifsh->pimpl->sys); if (ifsh->pimpl->nw > 0) { - ifsh_set_effective_well_params(WI, wdp, totmob, omega, ifsh); + ifsh_set_effective_well_params(WI, wdp, ifsh); hybsys_well_schur_comp_symm(ifsh->pimpl->nc, - ifsh->pimpl->cwpos, + ifsh->pimpl->cwell_pos, ifsh->pimpl->WI, ifsh->pimpl->sys, ifsh->pimpl->wsys); } - npp = ifsh_assemble_grid(bc, src, ifsh); + npp = ifsh_assemble_grid(bc, Binv, gpress, src, ifsh); if (ifsh->pimpl->nw > 0) { npp += ifsh_assemble_well(bc, wctrl, ifsh); @@ -634,7 +433,9 @@ ifsh_assemble(flowbc_t *bc, * substitution process, projected half-contact fluxes. */ /* ---------------------------------------------------------------------- */ void -ifsh_press_flux(grid_t *G, struct ifsh_data *h, +ifsh_press_flux(grid_t *G, + const double *Binv, const double *gpress, + struct fsh_data *h, double *cpress, double *fflux, double *wpress, double *wflux) /* ---------------------------------------------------------------------- */ @@ -645,8 +446,7 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h, hybsys_compute_press_flux(G->number_of_cells, G->cell_facepos, G->cell_faces, - h->pimpl->gpress, - h->pimpl->Binv, + gpress, Binv, h->pimpl->sys, h->x, cpress, h->pimpl->cflux, h->pimpl->work); @@ -655,8 +455,8 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h, assert ((wpress != NULL) && (wflux != NULL)); hybsys_compute_press_flux_well(G->number_of_cells, G->cell_facepos, G->number_of_faces, h->pimpl->nw, - h->pimpl->cwpos, h->pimpl->cwells, - h->pimpl->Binv, h->pimpl->WI, + h->pimpl->cwell_pos, h->pimpl->cwells, + Binv, h->pimpl->WI, h->pimpl->wdp, h->pimpl->sys, h->pimpl->wsys, h->x, cpress, h->pimpl->cflux, wpress, wflux, diff --git a/ifsh.h b/ifsh.h index cbefe5fd..8a51799e 100644 --- a/ifsh.h +++ b/ifsh.h @@ -23,8 +23,6 @@ #include "grid.h" #include "well.h" #include "flow_bc.h" -#include "sparse_sys.h" - #ifdef __cplusplus extern "C" { @@ -38,48 +36,19 @@ extern "C" { * and face fluxes is also included. */ +struct fsh_data; - -struct ifsh_impl; - -/** Contains the linear system for assembly, as well as internal data - * for the assembly routines. - */ -struct ifsh_data { - /* Let \f$n_i\f$ be the number of connections/faces of grid cell - * number \f$i\f$. Then max_ngconn = \f$\max_i n_i\f$ - */ - int max_ngconn; - /* With n_i as above, sum_ngconn2 = \f$\sum_i n_i^2\f$ */ - size_t sum_ngconn2; - - /* Linear system */ - struct CSRMatrix *A; /* Coefficient matrix */ - double *b; /* System RHS */ - double *x; /* Solution */ - - /* Private implementational details. */ - struct ifsh_impl *pimpl; -}; - - - -/** Constructs the ifsh_data object for a given grid and well - * pattern. +/** Constructs incompressible hybrid flow-solver data object for a + * given grid and well pattern. + * * @param G The grid * @param W The wells */ -struct ifsh_data * +struct fsh_data * ifsh_construct(grid_t *G, well_t *W); -/** Destroys the ifsh_data object */ -void -ifsh_destroy(struct ifsh_data *h); - - - /** Assembles the hybridized linear system for face pressures. * * This routine produces no output, other than changing the linear @@ -107,16 +76,14 @@ ifsh_destroy(struct ifsh_data *h); * be modified). Must already be constructed. */ void -ifsh_assemble(flowbc_t *bc, - double *src, - double *Binv, - double *gpress, - well_control_t *wctrl, - double *WI, - double *wdp, - double *totmob, /* \sum_i \lambda_i */ - double *omega, /* \sum_i \rho_i f_i */ - struct ifsh_data *h); +ifsh_assemble(flowbc_t *bc, + const double *src, + const double *Binv, + const double *gpress, + well_control_t *wctrl, + const double *WI, + const double *wdp, + struct fsh_data *h); /** Computes cell pressures, face fluxes, well pressures and well * fluxes from face pressures. @@ -131,7 +98,9 @@ ifsh_assemble(flowbc_t *bc, * @param wflux[out] \TODO */ void -ifsh_press_flux(grid_t *G, struct ifsh_data *h, +ifsh_press_flux(grid_t *G, + const double *Binv, const double *gpress, + struct fsh_data *h, double *cpress, double *fflux, double *wpress, double *wflux); From 4388dc68d6683604b16b4dbfdaa84c96bcfa5a19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Oct 2010 23:04:15 +0200 Subject: [PATCH 07/17] Add initial files for implementing compressible flows. --- Makefile.am | 30 ++-- fsh.c | 496 ++++++++++++++++++++++++++++++++++++++++++++++++++++ fsh.h | 71 ++++++++ 3 files changed, 583 insertions(+), 14 deletions(-) create mode 100644 fsh.c create mode 100644 fsh.h diff --git a/Makefile.am b/Makefile.am index 9fb8f72c..00ad0ae8 100644 --- a/Makefile.am +++ b/Makefile.am @@ -14,6 +14,7 @@ coarse_conn.h \ coarse_sys.h \ dfs.h \ flow_bc.h \ +fsh.h \ fsh_common.h \ grid.h \ hash_set.h \ @@ -29,20 +30,21 @@ GridAdapter.hpp \ HybridPressureSolver.hpp -libopmpressure_la_SOURCES = \ -coarse_conn.c \ -coarse_sys.c \ -dfs.c \ -flow_bc.c \ -fsh_common.c \ -hash_set.c \ -hybsys.c \ -hybsys_global.c \ -ifsh.c \ -ifsh_ms.c \ -mimetic.c \ -partition.c \ -sparse_sys.c \ +libopmpressure_la_SOURCES = \ +coarse_conn.c \ +coarse_sys.c \ +dfs.c \ +flow_bc.c \ +fsh.c \ +fsh_common.c \ +hash_set.c \ +hybsys.c \ +hybsys_global.c \ +ifsh.c \ +ifsh_ms.c \ +mimetic.c \ +partition.c \ +sparse_sys.c \ well.c diff --git a/fsh.c b/fsh.c new file mode 100644 index 00000000..d1f71f71 --- /dev/null +++ b/fsh.c @@ -0,0 +1,496 @@ +/* + 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 +#include +#include +#include +#include + +#include "fsh_common.h" +#include "ifsh.h" +#include "fsh_common_impl.h" +#include "hybsys.h" +#include "hybsys_global.h" + +#if defined MAX +#undef MAX +#endif +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + + +/* ---------------------------------------------------------------------- */ +static void +fsh_compute_table_sz(grid_t *G, well_t *W, int max_ngconn, + size_t *nnu, size_t *idata_sz, size_t *ddata_sz) +/* ---------------------------------------------------------------------- */ +{ + int nc, ngconn_tot; + + *nnu = G->number_of_faces; + + nc = G->number_of_cells; + ngconn_tot = G->cell_facepos[nc]; + + *idata_sz = nc + 1; /* gdof_pos */ + *idata_sz += ngconn_tot; /* gdof */ + *idata_sz += max_ngconn; /* iwork */ + + *ddata_sz = 2 * (*nnu); /* rhs + soln */ + *ddata_sz += ngconn_tot; /* cflux */ + *ddata_sz += max_ngconn; /* work */ + + if (W != NULL) { + *nnu += W->number_of_wells; + + /* cwell_pos */ + *idata_sz += nc + 1; + + /* cwells */ + *idata_sz += 2 * W->well_connpos[ W->number_of_wells ]; + + /* rhs + soln */ + *ddata_sz += 2 * W->number_of_wells; + + /* WI, wdp */ + *ddata_sz += 2 * W->well_connpos[ W->number_of_wells ]; + } +} + + +#if 0 +/* ---------------------------------------------------------------------- */ +static void +fsh_set_effective_well_params(const double *WI, + const double *wdp, + struct fsh_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, nc, i, perf; + int *cwpos, *cwells; + double *wsys_WI, *wsys_wdp; + + nc = h->pimpl->nc; + cwpos = h->pimpl->cwell_pos; + cwells = h->pimpl->cwells; + wsys_WI = h->pimpl->WI; + wsys_wdp = h->pimpl->wdp; + + for (c = i = 0; c < nc; c++) { + for (; i < cwpos[c + 1]; i++) { + perf = cwells[2*i + 1]; + + wsys_WI [i] = WI [perf]; + wsys_wdp[i] = wdp[perf]; + } + } +} +#endif + + +/* ---------------------------------------------------------------------- */ +static int +fsh_assemble_grid(flowbc_t *bc, + const double *Binv, + const double *gpress, + const double *src, + struct fsh_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, n, nc, p1, p2; + int npp; + int *pgconn, *gconn; + + nc = h->pimpl->nc; + pgconn = h->pimpl->gdof_pos; + gconn = h->pimpl->gdof; + + p1 = p2 = npp = 0; + for (c = 0; c < nc; c++) { + n = pgconn[c + 1] - pgconn[c]; + + hybsys_cellcontrib_unsymm(c, n, p1, p2, gpress, src, Binv, + h->pimpl->sys); + + npp += fsh_impose_bc(n, gconn + p1, bc, h->pimpl); + + hybsys_global_assemble_cell(n, gconn + p1, + h->pimpl->sys->S, + h->pimpl->sys->r, h->A, h->b); + + p1 += n; + p2 += n * n; + } + + return npp; +} + + +#if 0 +/* ---------------------------------------------------------------------- */ +static void +fsh_impose_well_control(int c, + flowbc_t *bc, + well_control_t *wctrl, + struct fsh_data *ifsh) +/* ---------------------------------------------------------------------- */ +{ + int ngconn, nwconn, i, w1, w2, wg, f; + int *pgconn, *gconn, *pwconn, *wconn; + + double bhp; + double *r, *r2w, *w2w; + + /* Enforce symmetric system */ + assert (ifsh->pimpl->wsys->r2w == ifsh->pimpl->wsys->w2r); + + pgconn = ifsh->pimpl->gdof_pos; + pwconn = ifsh->pimpl->cwell_pos; + + gconn = ifsh->pimpl->gdof + pgconn[c]; + wconn = ifsh->pimpl->cwells + 2*pwconn[c]; + + ngconn = pgconn[c + 1] - pgconn[c]; + nwconn = pwconn[c + 1] - pwconn[c]; + + r2w = ifsh->pimpl->wsys->r2w; + w2w = ifsh->pimpl->wsys->w2w; + r = ifsh->pimpl->wsys->r ; + + /* Adapt local system to prescribed boundary pressures (r->w) */ + for (i = 0; i < ngconn; i++) { + f = gconn[i]; + + if (bc->type[f] == PRESSURE) { + for (w1 = 0; w1 < nwconn; w1++) { + /* Eliminate prescribed (boundary) pressure value */ + r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->bcval[f]; + r2w[i + w1*ngconn] = 0.0; + } + + r[i] = 0.0; /* RHS value handled in *reservoir* asm */ + } + } + + /* Adapt local system to prescribed well (bottom-hole) pressures; + * w->r and w->w. */ + for (w1 = 0; w1 < nwconn; w1++) { + wg = wconn[2*w1 + 0]; + + if (wctrl->ctrl[wg] == BHP) { + bhp = wctrl->target[wg]; + + /* Well->reservoir */ + for (i = 0; i < ngconn; i++) { + assert ((bc->type[gconn[i]] != PRESSURE) || + !(fabs(r2w[i + w1*ngconn]) > 0.0)); + + r [i] -= r2w[i + w1*ngconn] * bhp; + r2w[i + w1*ngconn] = 0.0; + } + + /* Well->well */ + for (w2 = (w1 + 1) % nwconn; w2 != w1; w2 = (w2 + 1) % nwconn) { + r [ngconn + w2] -= w2w[w2 + w1*nwconn] * bhp; + w2w[w2 + w1*ngconn] = 0.0; + w2w[w1 + w2*ngconn] = 0.0; + } + + /* Assemble final well equation of the form S*p_bh = S*p_bh^0 */ + assert (fabs(w2w[w1 * (nwconn + 1)]) > 0.0); + + r[ngconn + w1] = w2w[w1 * (nwconn + 1)] * bhp; + } + } +} +#endif + + +#if 0 +/* ---------------------------------------------------------------------- */ +static int +fsh_assemble_well(flowbc_t *bc, + well_control_t *wctrl, + struct fsh_data *ifsh) +/* ---------------------------------------------------------------------- */ +{ + int npp; + int ngconn, nwconn, c, nc, w; + + int *pgconn, *gconn, *pwconn, *wconn; + + nc = ifsh->pimpl->nc; + + pgconn = ifsh->pimpl->gdof_pos; + gconn = ifsh->pimpl->gdof; + pwconn = ifsh->pimpl->cwell_pos; + wconn = ifsh->pimpl->cwells; + + for (c = 0; c < nc; c++) { + ngconn = pgconn[c + 1] - pgconn[c]; + nwconn = pwconn[c + 1] - pwconn[c]; + + if (nwconn > 0) { + hybsys_well_cellcontrib_symm(c, ngconn, pgconn[c], + pwconn, + ifsh->pimpl->WI, + ifsh->pimpl->wdp, + ifsh->pimpl->sys, + ifsh->pimpl->wsys); + + ifsh_impose_well_control(c, bc, wctrl, ifsh); + + hybsys_global_assemble_well_sym(ifsh->pimpl->nf, + ngconn, gconn + pgconn[c], + nwconn, wconn + 2*pwconn[c] + 0, + ifsh->pimpl->wsys->r2w, + ifsh->pimpl->wsys->w2w, + ifsh->pimpl->wsys->r, + ifsh->A, ifsh->b); + } + } + + npp = 0; + for (w = 0; w < ifsh->pimpl->nw; w++) { + if (wctrl->ctrl[w] == BHP) { + npp += 1; + } else if (wctrl->ctrl[w] == RATE) { + /* Impose total rate constraint. + * + * Note sign resulting from ->target[w] denoting + * *injection* flux. */ + ifsh->b[ifsh->pimpl->nf + w] -= - wctrl->target[w]; + } + } + + return npp; +} +#endif + + +/* ====================================================================== + * Public routines follow. + * ====================================================================== */ + + +/* ---------------------------------------------------------------------- */ +/* Allocate and define supporting structures for assembling the global + * system of linear equations to couple the grid (reservoir) + * connections represented by 'G' and, if present (i.e., non-NULL), + * the well connections represented by 'W'. */ +/* ---------------------------------------------------------------------- */ +struct fsh_data * +fsh_construct(grid_t *G, well_t *W) +/* ---------------------------------------------------------------------- */ +{ + int nc, ngconn_tot; + size_t idata_sz, ddata_sz, nnu; + struct fsh_data *new; + + assert (G != NULL); + + /* Allocate master structure, define system matrix sparsity */ + new = malloc(1 * sizeof *new); + if (new != NULL) { + new->A = hybsys_define_globconn(G, W); + new->pimpl = NULL; + + if (new->A == NULL) { + fsh_destroy(new); + new = NULL; + } + } + + + /* Allocate implementation structure */ + if (new != NULL) { + fsh_count_grid_dof(G, &new->max_ngconn, &new->sum_ngconn2); + + fsh_compute_table_sz(G, W, new->max_ngconn, + &nnu, &idata_sz, &ddata_sz); + + new->pimpl = fsh_impl_allocate_basic(idata_sz, ddata_sz); + + if (new->pimpl == NULL) { + fsh_destroy(new); + new = NULL; + } + } + + + /* Allocate Schur complement contributions. Symmetric system. */ + if (new != NULL) { + nc = G->number_of_cells; + ngconn_tot = G->cell_facepos[nc]; + + fsh_define_linsys_arrays(new); + fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn, + W, new->pimpl); + + new->pimpl->sys = hybsys_allocate_unsymm(new->max_ngconn, + nc, ngconn_tot); + + if (W != NULL) { + fsh_define_cell_wells(nc, W, new->pimpl); + + new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc, + new->pimpl->cwell_pos); + } + + if ((new->pimpl->sys == NULL) || + ((W != NULL) && (new->pimpl->wsys == NULL))) { + /* Failed to allocate ->sys or ->wsys (if W != NULL) */ + fsh_destroy(new); + new = NULL; + } + } + + + if (new != NULL) { + /* All allocations succeded. Fill metadata and return. */ + new->pimpl->nc = nc; + new->pimpl->nf = G->number_of_faces; + new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0; + + memcpy(new->pimpl->gdof_pos, + G->cell_facepos , + (nc + 1) * sizeof *new->pimpl->gdof_pos); + + memcpy(new->pimpl->gdof , + G->cell_faces , + ngconn_tot * sizeof *new->pimpl->gdof); + + hybsys_init(new->max_ngconn, new->pimpl->sys); + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +/* Assemble global system of linear equations + * + * ifsh->A * ifsh->x = ifsh->b + * + * from local inner product matrices Binv, gravity pressure gpress, + * boundary conditions bc, source terms src and fluid properties + * totmob and omega. */ +/* ---------------------------------------------------------------------- */ +void +fsh_assemble(flowbc_t *bc, + const double *src, + const double *Binv, + const double *Biv, + const double *P, + const double *gpress, + well_control_t *wctrl, + const double *WI, + const double *BivW, + const double *wdp, + struct fsh_data *h) +/* ---------------------------------------------------------------------- */ +{ + int npp; /* Number of prescribed pressure values */ + + hybsys_schur_comp_unsymm(h->pimpl->nc, + h->pimpl->gdof_pos, + Binv, Biv, P, h->pimpl->sys); + +#if 0 + if (ifsh->pimpl->nw > 0) { + ifsh_set_effective_well_params(WI, wdp, ifsh); + + hybsys_well_schur_comp_symm(ifsh->pimpl->nc, + ifsh->pimpl->cwell_pos, + ifsh->pimpl->WI, + ifsh->pimpl->sys, + ifsh->pimpl->wsys); + } +#endif + + npp = fsh_assemble_grid(bc, Binv, gpress, src, h); + +#if 0 + if (ifsh->pimpl->nw > 0) { + npp += ifsh_assemble_well(bc, wctrl, ifsh); + } +#endif + + if (npp == 0) { + h->A->sa[0] *= 2; /* Remove zero eigenvalue */ + } +} + + +/* ---------------------------------------------------------------------- */ +/* Compute cell pressures (cpress) and interface fluxes (fflux) from + * current solution of system of linear equations, h->x. Back + * substitution process, projected half-contact fluxes. */ +/* ---------------------------------------------------------------------- */ +void +ifsh_press_flux(grid_t *G, + const double *Binv, const double *gpress, + struct fsh_data *h, + double *cpress, double *fflux, + double *wpress, double *wflux) +/* ---------------------------------------------------------------------- */ +{ + int c, f, i; + double s; + + hybsys_compute_press_flux(G->number_of_cells, + G->cell_facepos, + G->cell_faces, + gpress, Binv, + h->pimpl->sys, + h->x, cpress, h->pimpl->cflux, + h->pimpl->work); + + if (h->pimpl->nw > 0) { + assert ((wpress != NULL) && (wflux != NULL)); + hybsys_compute_press_flux_well(G->number_of_cells, G->cell_facepos, + G->number_of_faces, h->pimpl->nw, + h->pimpl->cwell_pos, h->pimpl->cwells, + Binv, h->pimpl->WI, + h->pimpl->wdp, h->pimpl->sys, + h->pimpl->wsys, h->x, cpress, + h->pimpl->cflux, wpress, wflux, + h->pimpl->work); + } + + for (f = 0; f < G->number_of_faces; f++) { fflux[f] = 0.0; } + + i = 0; + for (c = 0; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0; + + fflux[f] += s * h->pimpl->cflux[i]; + } + } + + for (f = 0; f < G->number_of_faces; f++) { + i = (G->face_cells[2*f + 0] >= 0) + + (G->face_cells[2*f + 1] >= 0); + + fflux[f] /= i; + } +} diff --git a/fsh.h b/fsh.h new file mode 100644 index 00000000..4bc5b58b --- /dev/null +++ b/fsh.h @@ -0,0 +1,71 @@ +/* + 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_FSH_HEADER_INCLUDED +#define OPM_FHS_HEADER_INCLUDED + +#include "grid.h" +#include "well.h" +#include "flow_bc.h" + +#ifdef __cplusplus +extern "C" { +#endif + +struct fsh_data; + +/** Constructs incompressible hybrid flow-solver data object for a + * given grid and well pattern. + */ +struct fsh_data * +fsh_construct(grid_t *G, well_t *W); + + + +/** Assembles the hybridized linear system for face pressures. + */ +void +fsh_assemble(flowbc_t *bc, + const double *src, + const double *Binv, + const double *Biv, + const double *P, + const double *gpress, + well_control_t *wctrl, + const double *WI, + const double *BivW, + const double *wdp, + struct fsh_data *h); + +/** Computes cell pressures, face fluxes, well pressures and well + * fluxes from face pressures. + */ +void +fsh_press_flux(grid_t *G, + const double *Binv, const double *gpress, + struct fsh_data *h, + double *cpress, double *fflux, + double *wpress, double *wflux); + +#ifdef __cplusplus +} +#endif + + +#endif /* OPM_IFSH_HEADER_INCLUDED */ From 16378ac793c7cbe4356c4c7e6e7bdb46923b922c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Oct 2010 10:10:13 +0200 Subject: [PATCH 08/17] Janitorial work to re-enable building. Remove defunct well support inherited from ifsh.c. Rename ifsh_* functions to fsh_* functions. Update calls accordingly. --- fsh.c | 202 ++-------------------------------------------------------- 1 file changed, 6 insertions(+), 196 deletions(-) diff --git a/fsh.c b/fsh.c index d1f71f71..be35a0fb 100644 --- a/fsh.c +++ b/fsh.c @@ -25,7 +25,7 @@ #include #include "fsh_common.h" -#include "ifsh.h" +#include "fsh.h" #include "fsh_common_impl.h" #include "hybsys.h" #include "hybsys_global.h" @@ -75,36 +75,6 @@ fsh_compute_table_sz(grid_t *G, well_t *W, int max_ngconn, } -#if 0 -/* ---------------------------------------------------------------------- */ -static void -fsh_set_effective_well_params(const double *WI, - const double *wdp, - struct fsh_data *h) -/* ---------------------------------------------------------------------- */ -{ - int c, nc, i, perf; - int *cwpos, *cwells; - double *wsys_WI, *wsys_wdp; - - nc = h->pimpl->nc; - cwpos = h->pimpl->cwell_pos; - cwells = h->pimpl->cwells; - wsys_WI = h->pimpl->WI; - wsys_wdp = h->pimpl->wdp; - - for (c = i = 0; c < nc; c++) { - for (; i < cwpos[c + 1]; i++) { - perf = cwells[2*i + 1]; - - wsys_WI [i] = WI [perf]; - wsys_wdp[i] = wdp[perf]; - } - } -} -#endif - - /* ---------------------------------------------------------------------- */ static int fsh_assemble_grid(flowbc_t *bc, @@ -143,148 +113,6 @@ fsh_assemble_grid(flowbc_t *bc, } -#if 0 -/* ---------------------------------------------------------------------- */ -static void -fsh_impose_well_control(int c, - flowbc_t *bc, - well_control_t *wctrl, - struct fsh_data *ifsh) -/* ---------------------------------------------------------------------- */ -{ - int ngconn, nwconn, i, w1, w2, wg, f; - int *pgconn, *gconn, *pwconn, *wconn; - - double bhp; - double *r, *r2w, *w2w; - - /* Enforce symmetric system */ - assert (ifsh->pimpl->wsys->r2w == ifsh->pimpl->wsys->w2r); - - pgconn = ifsh->pimpl->gdof_pos; - pwconn = ifsh->pimpl->cwell_pos; - - gconn = ifsh->pimpl->gdof + pgconn[c]; - wconn = ifsh->pimpl->cwells + 2*pwconn[c]; - - ngconn = pgconn[c + 1] - pgconn[c]; - nwconn = pwconn[c + 1] - pwconn[c]; - - r2w = ifsh->pimpl->wsys->r2w; - w2w = ifsh->pimpl->wsys->w2w; - r = ifsh->pimpl->wsys->r ; - - /* Adapt local system to prescribed boundary pressures (r->w) */ - for (i = 0; i < ngconn; i++) { - f = gconn[i]; - - if (bc->type[f] == PRESSURE) { - for (w1 = 0; w1 < nwconn; w1++) { - /* Eliminate prescribed (boundary) pressure value */ - r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->bcval[f]; - r2w[i + w1*ngconn] = 0.0; - } - - r[i] = 0.0; /* RHS value handled in *reservoir* asm */ - } - } - - /* Adapt local system to prescribed well (bottom-hole) pressures; - * w->r and w->w. */ - for (w1 = 0; w1 < nwconn; w1++) { - wg = wconn[2*w1 + 0]; - - if (wctrl->ctrl[wg] == BHP) { - bhp = wctrl->target[wg]; - - /* Well->reservoir */ - for (i = 0; i < ngconn; i++) { - assert ((bc->type[gconn[i]] != PRESSURE) || - !(fabs(r2w[i + w1*ngconn]) > 0.0)); - - r [i] -= r2w[i + w1*ngconn] * bhp; - r2w[i + w1*ngconn] = 0.0; - } - - /* Well->well */ - for (w2 = (w1 + 1) % nwconn; w2 != w1; w2 = (w2 + 1) % nwconn) { - r [ngconn + w2] -= w2w[w2 + w1*nwconn] * bhp; - w2w[w2 + w1*ngconn] = 0.0; - w2w[w1 + w2*ngconn] = 0.0; - } - - /* Assemble final well equation of the form S*p_bh = S*p_bh^0 */ - assert (fabs(w2w[w1 * (nwconn + 1)]) > 0.0); - - r[ngconn + w1] = w2w[w1 * (nwconn + 1)] * bhp; - } - } -} -#endif - - -#if 0 -/* ---------------------------------------------------------------------- */ -static int -fsh_assemble_well(flowbc_t *bc, - well_control_t *wctrl, - struct fsh_data *ifsh) -/* ---------------------------------------------------------------------- */ -{ - int npp; - int ngconn, nwconn, c, nc, w; - - int *pgconn, *gconn, *pwconn, *wconn; - - nc = ifsh->pimpl->nc; - - pgconn = ifsh->pimpl->gdof_pos; - gconn = ifsh->pimpl->gdof; - pwconn = ifsh->pimpl->cwell_pos; - wconn = ifsh->pimpl->cwells; - - for (c = 0; c < nc; c++) { - ngconn = pgconn[c + 1] - pgconn[c]; - nwconn = pwconn[c + 1] - pwconn[c]; - - if (nwconn > 0) { - hybsys_well_cellcontrib_symm(c, ngconn, pgconn[c], - pwconn, - ifsh->pimpl->WI, - ifsh->pimpl->wdp, - ifsh->pimpl->sys, - ifsh->pimpl->wsys); - - ifsh_impose_well_control(c, bc, wctrl, ifsh); - - hybsys_global_assemble_well_sym(ifsh->pimpl->nf, - ngconn, gconn + pgconn[c], - nwconn, wconn + 2*pwconn[c] + 0, - ifsh->pimpl->wsys->r2w, - ifsh->pimpl->wsys->w2w, - ifsh->pimpl->wsys->r, - ifsh->A, ifsh->b); - } - } - - npp = 0; - for (w = 0; w < ifsh->pimpl->nw; w++) { - if (wctrl->ctrl[w] == BHP) { - npp += 1; - } else if (wctrl->ctrl[w] == RATE) { - /* Impose total rate constraint. - * - * Note sign resulting from ->target[w] denoting - * *injection* flux. */ - ifsh->b[ifsh->pimpl->nf + w] -= - wctrl->target[w]; - } - } - - return npp; -} -#endif - - /* ====================================================================== * Public routines follow. * ====================================================================== */ @@ -413,26 +241,8 @@ fsh_assemble(flowbc_t *bc, h->pimpl->gdof_pos, Binv, Biv, P, h->pimpl->sys); -#if 0 - if (ifsh->pimpl->nw > 0) { - ifsh_set_effective_well_params(WI, wdp, ifsh); - - hybsys_well_schur_comp_symm(ifsh->pimpl->nc, - ifsh->pimpl->cwell_pos, - ifsh->pimpl->WI, - ifsh->pimpl->sys, - ifsh->pimpl->wsys); - } -#endif - npp = fsh_assemble_grid(bc, Binv, gpress, src, h); -#if 0 - if (ifsh->pimpl->nw > 0) { - npp += ifsh_assemble_well(bc, wctrl, ifsh); - } -#endif - if (npp == 0) { h->A->sa[0] *= 2; /* Remove zero eigenvalue */ } @@ -445,11 +255,11 @@ fsh_assemble(flowbc_t *bc, * substitution process, projected half-contact fluxes. */ /* ---------------------------------------------------------------------- */ void -ifsh_press_flux(grid_t *G, - const double *Binv, const double *gpress, - struct fsh_data *h, - double *cpress, double *fflux, - double *wpress, double *wflux) +fsh_press_flux(grid_t *G, + const double *Binv, const double *gpress, + struct fsh_data *h, + double *cpress, double *fflux, + double *wpress, double *wflux) /* ---------------------------------------------------------------------- */ { int c, f, i; From a88616529086b7a45e6d6ac02a49897917e2283e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Oct 2010 10:11:36 +0200 Subject: [PATCH 09/17] Correct spelling error in comment. --- fsh.c | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/fsh.c b/fsh.c index be35a0fb..ec6f3e07 100644 --- a/fsh.c +++ b/fsh.c @@ -215,11 +215,8 @@ fsh_construct(grid_t *G, well_t *W) /* ---------------------------------------------------------------------- */ /* Assemble global system of linear equations * - * ifsh->A * ifsh->x = ifsh->b - * - * from local inner product matrices Binv, gravity pressure gpress, - * boundary conditions bc, source terms src and fluid properties - * totmob and omega. */ + * fsh->A * fsh->x = fsh->b + */ /* ---------------------------------------------------------------------- */ void fsh_assemble(flowbc_t *bc, From b193e461cac687b37f0947c4d092b51fc69a1fe7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Oct 2010 19:14:07 +0200 Subject: [PATCH 10/17] Add utilities for computing derived fluid quantities. These are (currently) entirely geared towards pressure solvers. --- Makefile.am | 32 ++++++++++---------- compr_quant.c | 81 +++++++++++++++++++++++++++++++++++++++++++++++++++ compr_quant.h | 46 +++++++++++++++++++++++++++++ 3 files changed, 144 insertions(+), 15 deletions(-) create mode 100644 compr_quant.c create mode 100644 compr_quant.h diff --git a/Makefile.am b/Makefile.am index 00ad0ae8..8655e018 100644 --- a/Makefile.am +++ b/Makefile.am @@ -12,6 +12,7 @@ libopmpressure_HEADERS = \ blas_lapack.h \ coarse_conn.h \ coarse_sys.h \ +compr_quant.h \ dfs.h \ flow_bc.h \ fsh.h \ @@ -30,21 +31,22 @@ GridAdapter.hpp \ HybridPressureSolver.hpp -libopmpressure_la_SOURCES = \ -coarse_conn.c \ -coarse_sys.c \ -dfs.c \ -flow_bc.c \ -fsh.c \ -fsh_common.c \ -hash_set.c \ -hybsys.c \ -hybsys_global.c \ -ifsh.c \ -ifsh_ms.c \ -mimetic.c \ -partition.c \ -sparse_sys.c \ +libopmpressure_la_SOURCES = \ +coarse_conn.c \ +coarse_sys.c \ +compr_quant.c \ +dfs.c \ +flow_bc.c \ +fsh.c \ +fsh_common.c \ +hash_set.c \ +hybsys.c \ +hybsys_global.c \ +ifsh.c \ +ifsh_ms.c \ +mimetic.c \ +partition.c \ +sparse_sys.c \ well.c diff --git a/compr_quant.c b/compr_quant.c new file mode 100644 index 00000000..296168c7 --- /dev/null +++ b/compr_quant.c @@ -0,0 +1,81 @@ +/* + 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 "compr_quant.h" + + +/* ---------------------------------------------------------------------- */ +/* 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/compr_quant.h b/compr_quant.h new file mode 100644 index 00000000..fe269efd --- /dev/null +++ b/compr_quant.h @@ -0,0 +1,46 @@ +/* + 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" + +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); + +#endif /* OPM_COMPR_QUANT_HEADER_INCLUDED */ From 3715c2b4039ba0fd2df19d504283f7acb40bd714 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Oct 2010 19:17:58 +0200 Subject: [PATCH 11/17] Trivial, janitorial change. Correct a comment and break a long line. --- fsh.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/fsh.c b/fsh.c index ec6f3e07..428aa02b 100644 --- a/fsh.c +++ b/fsh.c @@ -163,7 +163,7 @@ fsh_construct(grid_t *G, well_t *W) } - /* Allocate Schur complement contributions. Symmetric system. */ + /* Allocate Schur complement contributions. Unsymmetric system. */ if (new != NULL) { nc = G->number_of_cells; ngconn_tot = G->cell_facepos[nc]; @@ -178,8 +178,9 @@ fsh_construct(grid_t *G, well_t *W) if (W != NULL) { fsh_define_cell_wells(nc, W, new->pimpl); - new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc, - new->pimpl->cwell_pos); + new->pimpl->wsys = + hybsys_well_allocate_unsymm(new->max_ngconn, nc, + new->pimpl->cwell_pos); } if ((new->pimpl->sys == NULL) || From 98afc235bc306f113d02b5fd5f5d2a18b9e23a8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 21 Oct 2010 14:13:17 +0200 Subject: [PATCH 12/17] Change 'P' term sign in accordance with MRST impl. --- compr_quant.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compr_quant.c b/compr_quant.c index 296168c7..30edf955 100644 --- a/compr_quant.c +++ b/compr_quant.c @@ -58,7 +58,7 @@ compr_accum_term(size_t nc, size_t c; for (c = 0; c < nc; c++) { - P[c] = totcompr[c] * porevol[c] / dt; + P[c] = -totcompr[c] * porevol[c] / dt; } } @@ -76,6 +76,6 @@ compr_src_add_press_accum(size_t nc, size_t c; for (c = 0; c < nc; c++) { - src[c] += P[c] * p0[c]; + src[c] -= P[c] * p0[c]; } } From 58591b7503a3ab48b9dcffc3d38652125d582384 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 24 Oct 2010 15:26:21 +0200 Subject: [PATCH 13/17] Add preliminary TPFA transmissibility calculator. --- Makefile.am | 2 ++ trans_tpfa.c | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++ trans_tpfa.h | 31 ++++++++++++++++++++ 3 files changed, 115 insertions(+) create mode 100644 trans_tpfa.c create mode 100644 trans_tpfa.h diff --git a/Makefile.am b/Makefile.am index 8655e018..160e120b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -26,6 +26,7 @@ ifsh_ms.h \ mimetic.h \ partition.h \ sparse_sys.h \ +trans_tpfa.h \ well.h \ GridAdapter.hpp \ HybridPressureSolver.hpp @@ -47,6 +48,7 @@ ifsh_ms.c \ mimetic.c \ partition.c \ sparse_sys.c \ +trans_tpfa.c \ well.c diff --git a/trans_tpfa.c b/trans_tpfa.c new file mode 100644 index 00000000..611d0387 --- /dev/null +++ b/trans_tpfa.c @@ -0,0 +1,82 @@ +#include +#include + +#include "blas_lapack.h" +#include "trans_tpfa.h" + +/* ---------------------------------------------------------------------- */ +/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */ +/* ---------------------------------------------------------------------- */ +void +tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans) +/* ---------------------------------------------------------------------- */ +{ + int c, d, f, i, j; + double s, dist, denom; + + double Kn[3]; + double *cc, *fc, *n; + double *K; + + MAT_SIZE_T nrows, ncols, ldA, incx, incy; + double a1, a2; + + d = G->dimensions; + + nrows = ncols = ldA = d; + incx = incy = 1 ; + a1 = 1.0; a2 = 0.0 ; + + for (c = i = 0; c < G->number_of_cells; c++) { + K = perm + (c * d * d); + cc = G->cell_centroids + (c * d); + + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0; + + n = G->face_normals + (f * d); + fc = G->face_centroids + (f * d); + + dgemv_("No Transpose", &nrows, &ncols, + &a1, K, &ldA, n, &incx, &a2, &Kn[0], &incy); + + htrans[i] = denom = 0.0; + for (j = 0; j < d; j++) { + dist = fc[j] - cc[j]; + + htrans[i] += s * dist * Kn[j]; + denom += dist * dist; + } + + assert (denom > 0); + htrans[i] /= denom; + htrans[i] = fabs(htrans[i]); + } + } +} + + +/* ---------------------------------------------------------------------- */ +void +tpfa_trans_compute(grid_t *G, const double *htrans, double *trans) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f; + + for (f = 0; f < G->number_of_faces; f++) { + trans[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]; + + trans[f] += 1.0 / htrans[i]; + } + } + + for (f = 0; f < G->number_of_faces; f++) { + trans[f] = 1.0 / trans[f]; + } +} diff --git a/trans_tpfa.h b/trans_tpfa.h new file mode 100644 index 00000000..ea08bfd9 --- /dev/null +++ b/trans_tpfa.h @@ -0,0 +1,31 @@ +/* + 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_TRANS_TPFA_HEADER_INCLUDED +#define OPM_TRANS_TPFA_HEADER_INCLUDED + +#include "grid.h" + +void +tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans); + +void +tpfa_trans_compute(grid_t *G, const double *htrans, double *trans); + +#endif /* OPM_TRANS_TPFA_HEADER_INCLUDED */ From afc77d25d147646914effe3034e628adcf4b7089 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 24 Oct 2010 21:20:15 +0200 Subject: [PATCH 14/17] Add (incomp) pressure system assembler based on TPFA. --- Makefile.am | 2 + ifs_tpfa.c | 227 ++++++++++++++++++++++++++++++++++++++++++++++++++++ ifs_tpfa.h | 57 +++++++++++++ 3 files changed, 286 insertions(+) create mode 100644 ifs_tpfa.c create mode 100644 ifs_tpfa.h diff --git a/Makefile.am b/Makefile.am index 160e120b..f2543169 100644 --- a/Makefile.am +++ b/Makefile.am @@ -21,6 +21,7 @@ grid.h \ hash_set.h \ hybsys_global.h \ hybsys.h \ +ifs_tpfa.h \ ifsh.h \ ifsh_ms.h \ mimetic.h \ @@ -43,6 +44,7 @@ fsh_common.c \ hash_set.c \ hybsys.c \ hybsys_global.c \ +ifs_tpfa.c \ ifsh.c \ ifsh_ms.c \ mimetic.c \ diff --git a/ifs_tpfa.c b/ifs_tpfa.c new file mode 100644 index 00000000..419b1597 --- /dev/null +++ b/ifs_tpfa.c @@ -0,0 +1,227 @@ +#include +#include +#include +#include + +#include "ifs_tpfa.h" +#include "sparse_sys.h" + + +struct ifs_tpfa_impl { + /* Linear storage */ + double *ddata; +}; + + +/* ---------------------------------------------------------------------- */ +static void +impl_deallocate(struct ifs_tpfa_impl *pimpl) +/* ---------------------------------------------------------------------- */ +{ + if (pimpl != NULL) { + free(pimpl->ddata); + } + + free(pimpl); +} + + +/* ---------------------------------------------------------------------- */ +static struct ifs_tpfa_impl * +impl_allocate(grid_t *G) +/* ---------------------------------------------------------------------- */ +{ + struct ifs_tpfa_impl *new; + + size_t ddata_sz; + + ddata_sz = 2 * G->number_of_cells; /* b, x */ + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->ddata = malloc(ddata_sz * sizeof *new->ddata); + + if (new->ddata == NULL) { + impl_deallocate(new); + new = NULL; + } + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +static struct CSRMatrix * +ifs_tpfa_construct_matrix(grid_t *G) +/* ---------------------------------------------------------------------- */ +{ + int f, c1, c2; + size_t nnz; + + struct CSRMatrix *A; + + A = csrmatrix_new_count_nnz(G->number_of_cells); + + if (A != NULL) { + /* Self connections */ + for (c1 = 0; c1 < G->number_of_cells; c1++) { + A->ia[ c1 + 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; + } + } + + nnz = csrmatrix_new_elms_pushback(A); + if (nnz == 0) { + csrmatrix_delete(A); + A = NULL; + } + } + + if (A != NULL) { + /* Fill self connections */ + for (c1 = 0; c1 < G->number_of_cells; c1++) { + A->ja[ A->ia[ c1 + 1 ] ++ ] = c1; + } + + /* 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; + } + } + + assert ((size_t) A->ia[ G->number_of_cells ] == nnz); + + /* Guarantee sorted rows */ + csrmatrix_sortrows(A); + } + + return A; +} + + + +/* ====================================================================== + * Public interface below separator. + * ====================================================================== */ + +/* ---------------------------------------------------------------------- */ +struct ifs_tpfa_data * +ifs_tpfa_construct(grid_t *G) +/* ---------------------------------------------------------------------- */ +{ + struct ifs_tpfa_data *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->pimpl = impl_allocate(G); + new->A = ifs_tpfa_construct_matrix(G); + + if ((new->pimpl == NULL) || (new->A == NULL)) { + ifs_tpfa_destroy(new); + new = NULL; + } + } + + if (new != NULL) { + new->b = new->pimpl->ddata; + new->x = new->b + new->A->m; + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +void +ifs_tpfa_assemble(grid_t *G, + const double *trans, + const double *src, + struct ifs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2, c, i, f, j1, j2; + + csrmatrix_zero( h->A); + vector_zero (h->A->m, h->b); + + for (c = i = 0; c < G->number_of_cells; c++) { + j1 = csrmatrix_elm_index(c, c, h->A); + + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; + + c2 = (c1 == c) ? c2 : c1; + + if (c2 >= 0) { + j2 = csrmatrix_elm_index(c, c2, h->A); + + h->A->sa[j1] += trans[f]; + h->A->sa[j2] -= trans[f]; + } + } + + h->b[c] += src[c]; + } +} + + +/* ---------------------------------------------------------------------- */ +void +ifs_tpfa_press_flux(grid_t *G, + const double *trans, + const double *src, + struct ifs_tpfa_data *h, + double *cpress, + double *fflux) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2, f; + + /* Assign cell pressure directly from solution vector */ + memcpy(cpress, h->x, G->number_of_cells * sizeof *cpress); + + 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)) { + fflux[f] = trans[f] * (cpress[c1] - cpress[c2]); + } else { + fflux[f] = 0.0; + } + } +} + + +/* ---------------------------------------------------------------------- */ +void +ifs_tpfa_destroy(struct ifs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + if (h != NULL) { + csrmatrix_delete(h->A); + impl_deallocate (h->pimpl); + } + + free(h); +} diff --git a/ifs_tpfa.h b/ifs_tpfa.h new file mode 100644 index 00000000..20933c1c --- /dev/null +++ b/ifs_tpfa.h @@ -0,0 +1,57 @@ +/* + 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_IFS_TPFA_HEADER_INCLUDED +#define OPM_IFS_TPFA_HEADER_INCLUDED + +#include "grid.h" + +struct ifs_tpfa_impl; +struct CSRMatrix; + +struct ifs_tpfa_data { + struct CSRMatrix *A; + double *b; + double *x; + + struct ifs_tpfa_impl *pimpl; +}; + + +struct ifs_tpfa_data * +ifs_tpfa_construct(grid_t *G); + +void +ifs_tpfa_assemble(grid_t *G, + const double *trans, + const double *src, + struct ifs_tpfa_data *h); + +void +ifs_tpfa_press_flux(grid_t *G, + const double *trans, + const double *src, + struct ifs_tpfa_data *h, + double *cpress, + double *fflux); + +void +ifs_tpfa_destroy(struct ifs_tpfa_data *h); + +#endif /* OPM_IFS_TPFA_HEADER_INCLUDED */ From 1965507cad42dbd6bf3f50b728dc48d27737dc7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 24 Oct 2010 21:22:29 +0200 Subject: [PATCH 15/17] The ifs_tpfa_press_flux() function does not need 'src'. --- ifs_tpfa.c | 1 - ifs_tpfa.h | 1 - 2 files changed, 2 deletions(-) diff --git a/ifs_tpfa.c b/ifs_tpfa.c index 419b1597..5a633fca 100644 --- a/ifs_tpfa.c +++ b/ifs_tpfa.c @@ -189,7 +189,6 @@ ifs_tpfa_assemble(grid_t *G, void ifs_tpfa_press_flux(grid_t *G, const double *trans, - const double *src, struct ifs_tpfa_data *h, double *cpress, double *fflux) diff --git a/ifs_tpfa.h b/ifs_tpfa.h index 20933c1c..433e9bcf 100644 --- a/ifs_tpfa.h +++ b/ifs_tpfa.h @@ -46,7 +46,6 @@ ifs_tpfa_assemble(grid_t *G, void ifs_tpfa_press_flux(grid_t *G, const double *trans, - const double *src, struct ifs_tpfa_data *h, double *cpress, double *fflux); From f486eb01fe21d056d64639e927c039005f1b9c0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 25 Oct 2010 13:56:38 +0200 Subject: [PATCH 16/17] Add traditional trans mobility weighting. --- trans_tpfa.c | 28 ++++++++++++++++++++++++++++ trans_tpfa.h | 6 ++++++ 2 files changed, 34 insertions(+) diff --git a/trans_tpfa.c b/trans_tpfa.c index 611d0387..86cfa75c 100644 --- a/trans_tpfa.c +++ b/trans_tpfa.c @@ -80,3 +80,31 @@ tpfa_trans_compute(grid_t *G, const double *htrans, double *trans) trans[f] = 1.0 / trans[f]; } } + + +/* ---------------------------------------------------------------------- */ +void +tpfa_eff_trans_compute(grid_t *G, + const double *totmob, + const double *htrans, + double *trans) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f; + + for (f = 0; f < G->number_of_faces; f++) { + trans[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]; + + trans[f] += 1.0 / (totmob[c] * htrans[i]); + } + } + + for (f = 0; f < G->number_of_faces; f++) { + trans[f] = 1.0 / trans[f]; + } +} diff --git a/trans_tpfa.h b/trans_tpfa.h index ea08bfd9..a5575dda 100644 --- a/trans_tpfa.h +++ b/trans_tpfa.h @@ -28,4 +28,10 @@ tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans); void tpfa_trans_compute(grid_t *G, const double *htrans, double *trans); +void +tpfa_eff_trans_compute(grid_t *G, + const double *totmob, + const double *htrans, + double *trans); + #endif /* OPM_TRANS_TPFA_HEADER_INCLUDED */ From db0517f794879f5c348058d2cf7ba2a612977f0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 25 Oct 2010 20:28:34 +0200 Subject: [PATCH 17/17] Remove zero eigenval. Prepare for adding gravity. --- ifs_tpfa.c | 13 +++++++++++-- ifs_tpfa.h | 1 + 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/ifs_tpfa.c b/ifs_tpfa.c index 5a633fca..61608840 100644 --- a/ifs_tpfa.c +++ b/ifs_tpfa.c @@ -8,6 +8,8 @@ struct ifs_tpfa_impl { + double *fgrav; /* Accumulated grav contrib/face */ + /* Linear storage */ double *ddata; }; @@ -36,6 +38,7 @@ impl_allocate(grid_t *G) size_t ddata_sz; ddata_sz = 2 * G->number_of_cells; /* b, x */ + ddata_sz += 1 * G->number_of_faces; /* fgrav */ new = malloc(1 * sizeof *new); @@ -142,6 +145,8 @@ ifs_tpfa_construct(grid_t *G) if (new != NULL) { new->b = new->pimpl->ddata; new->x = new->b + new->A->m; + + new->pimpl->fgrav = new->x + new->A->m; } return new; @@ -153,13 +158,15 @@ void ifs_tpfa_assemble(grid_t *G, const double *trans, const double *src, + const double *gpress, struct ifs_tpfa_data *h) /* ---------------------------------------------------------------------- */ { int c1, c2, c, i, f, j1, j2; - csrmatrix_zero( h->A); - vector_zero (h->A->m, h->b); + csrmatrix_zero( h->A); + vector_zero (h->A->m, h->b); + vector_zero (G->number_of_faces, h->pimpl->fgrav); for (c = i = 0; c < G->number_of_cells; c++) { j1 = csrmatrix_elm_index(c, c, h->A); @@ -182,6 +189,8 @@ ifs_tpfa_assemble(grid_t *G, h->b[c] += src[c]; } + + h->A->sa[0] *= 2; } diff --git a/ifs_tpfa.h b/ifs_tpfa.h index 433e9bcf..4d4b08ca 100644 --- a/ifs_tpfa.h +++ b/ifs_tpfa.h @@ -41,6 +41,7 @@ void ifs_tpfa_assemble(grid_t *G, const double *trans, const double *src, + const double *gpress, struct ifs_tpfa_data *h); void