Add (incomp) pressure system assembler based on TPFA.

ifs_tpfa.c

@@ -0,0 +1,227 @@
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#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);
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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 */