diff --git a/Makefile.am b/Makefile.am index 40748115..f2543169 100644 --- a/Makefile.am +++ b/Makefile.am @@ -12,17 +12,22 @@ libopmpressure_HEADERS = \ blas_lapack.h \ coarse_conn.h \ coarse_sys.h \ +compr_quant.h \ dfs.h \ flow_bc.h \ +fsh.h \ +fsh_common.h \ grid.h \ hash_set.h \ hybsys_global.h \ hybsys.h \ +ifs_tpfa.h \ ifsh.h \ ifsh_ms.h \ mimetic.h \ partition.h \ sparse_sys.h \ +trans_tpfa.h \ well.h \ GridAdapter.hpp \ HybridPressureSolver.hpp @@ -31,16 +36,21 @@ HybridPressureSolver.hpp 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 \ +ifs_tpfa.c \ ifsh.c \ ifsh_ms.c \ mimetic.c \ partition.c \ sparse_sys.c \ +trans_tpfa.c \ well.c 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.c b/coarse_sys.c index 793e632d..780cffd8 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++; } } @@ -1053,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/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/compr_quant.c b/compr_quant.c new file mode 100644 index 00000000..30edf955 --- /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 */ 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/fsh.c b/fsh.c new file mode 100644 index 00000000..428aa02b --- /dev/null +++ b/fsh.c @@ -0,0 +1,304 @@ +/* + 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 "fsh.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 ]; + } +} + + +/* ---------------------------------------------------------------------- */ +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; +} + + +/* ====================================================================== + * 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. Unsymmetric 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_unsymm(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 + * + * fsh->A * fsh->x = fsh->b + */ +/* ---------------------------------------------------------------------- */ +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); + + npp = fsh_assemble_grid(bc, Binv, gpress, src, h); + + 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 +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; + 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 */ 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/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/ifs_tpfa.c b/ifs_tpfa.c new file mode 100644 index 00000000..61608840 --- /dev/null +++ b/ifs_tpfa.c @@ -0,0 +1,235 @@ +#include +#include +#include +#include + +#include "ifs_tpfa.h" +#include "sparse_sys.h" + + +struct ifs_tpfa_impl { + double *fgrav; /* Accumulated grav contrib/face */ + + /* 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 */ + ddata_sz += 1 * G->number_of_faces; /* fgrav */ + + 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; + + new->pimpl->fgrav = new->x + new->A->m; + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +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); + 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); + + 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]; + } + + h->A->sa[0] *= 2; +} + + +/* ---------------------------------------------------------------------- */ +void +ifs_tpfa_press_flux(grid_t *G, + const double *trans, + 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..4d4b08ca --- /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, + const double *gpress, + struct ifs_tpfa_data *h); + +void +ifs_tpfa_press_flux(grid_t *G, + const double *trans, + struct ifs_tpfa_data *h, + double *cpress, + double *fflux); + +void +ifs_tpfa_destroy(struct ifs_tpfa_data *h); + +#endif /* OPM_IFS_TPFA_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 2f0a0dfe..8a51799e 100644 --- a/ifsh.h +++ b/ifsh.h @@ -17,14 +17,12 @@ 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" #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); @@ -141,4 +110,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.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/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.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 1b73aa1c..82e347cf 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" { @@ -81,8 +81,24 @@ 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 -#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.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 e9b9d75d..aadaedf5 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 @@ -72,8 +72,14 @@ 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 -#endif /* SPARSE_SYS_H_INCLUDED */ +#endif /* OPM_SPARSE_SYS_HEADER_INCLUDED */ diff --git a/trans_tpfa.c b/trans_tpfa.c new file mode 100644 index 00000000..86cfa75c --- /dev/null +++ b/trans_tpfa.c @@ -0,0 +1,110 @@ +#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]; + } +} + + +/* ---------------------------------------------------------------------- */ +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 new file mode 100644 index 00000000..a5575dda --- /dev/null +++ b/trans_tpfa.h @@ -0,0 +1,37 @@ +/* + 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); + +void +tpfa_eff_trans_compute(grid_t *G, + const double *totmob, + const double *htrans, + double *trans); + +#endif /* OPM_TRANS_TPFA_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 */