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] 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 */