Add files for implementing compressible TPFA assembly.

src/cfs_tpfa.c

@@ -0,0 +1,243 @@
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cfs_tpfa.h"
+#include "sparse_sys.h"
+
+
+struct cfs_tpfa_impl {
+    double *fgrav;              /* Accumulated grav contrib/face */
+
+    /* Linear storage */
+    double *ddata;
+};
+
+
+/* ---------------------------------------------------------------------- */
+static void
+impl_deallocate(struct cfs_tpfa_impl *pimpl)
+/* ---------------------------------------------------------------------- */
+{
+    if (pimpl != NULL) {
+        free(pimpl->ddata);
+    }
+
+    free(pimpl);
+}
+
+
+/* ---------------------------------------------------------------------- */
+static struct cfs_tpfa_impl *
+impl_allocate(grid_t *G)
+/* ---------------------------------------------------------------------- */
+{
+    struct cfs_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 *
+cfs_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;
+            }
+        }
+
+        /* The tpfa matrix is square */
+        A->n = A->m;
+
+        assert ((size_t) A->ia[ G->number_of_cells ] == nnz);
+
+        /* Guarantee sorted rows */
+        csrmatrix_sortrows(A);
+    }
+
+    return A;
+}
+
+
+
+/* ======================================================================
+ * Public interface below separator.
+ * ====================================================================== */
+
+/* ---------------------------------------------------------------------- */
+struct cfs_tpfa_data *
+cfs_tpfa_construct(grid_t *G)
+/* ---------------------------------------------------------------------- */
+{
+    struct cfs_tpfa_data *new;
+
+    new = malloc(1 * sizeof *new);
+
+    if (new != NULL) {
+        new->pimpl = impl_allocate(G);
+        new->A     = cfs_tpfa_construct_matrix(G);
+
+        if ((new->pimpl == NULL) || (new->A == NULL)) {
+            cfs_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
+cfs_tpfa_assemble(grid_t               *G,
+                  const double         *ctrans,
+                  const double         *P,
+                  const double         *src,
+                  struct cfs_tpfa_data *h)
+/* ---------------------------------------------------------------------- */
+{
+    int c1, c2, c, i, f, j1, j2;
+
+    double s;
+
+    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];
+
+            s  = 2.0*(c1 == c) - 1.0;
+            c2 = (c1 == c) ? c2 : c1;
+
+            if (c2 >= 0) {
+                j2 = csrmatrix_elm_index(c, c2, h->A);
+
+                h->A->sa[j1] += ctrans[i];
+                h->A->sa[j2] -= ctrans[i];
+            }
+        }
+
+        h->b[c] += src[c];
+
+        /* Compressible accumulation term */
+        h->A->sa[j1] += P[c];
+    }
+
+    h->A->sa[0] *= 2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+cfs_tpfa_press_flux(grid_t               *G,
+                    const double         *trans,
+                    struct cfs_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
+cfs_tpfa_destroy(struct cfs_tpfa_data *h)
+/* ---------------------------------------------------------------------- */
+{
+    if (h != NULL) {
+        csrmatrix_delete(h->A);
+        impl_deallocate (h->pimpl);
+    }
+
+    free(h);
+}