Add demonstration for bare-bones reordering implementation.

.hgignore

@@ -33,6 +33,7 @@ missing
 *_test
 examples/scaneclipsedeck
 examples/spu_2p
+examples/reorder-qfs
 tests/test_cfs_tpfa
 tests/test_jacsys
 tests/test_readvector

examples/Makefile.am

@@ -10,6 +10,6 @@ noinst_PROGRAMS =				\
 scaneclipsedeck
 
 if UMFPACK
-noinst_PROGRAMS += spu_2p
+noinst_PROGRAMS += spu_2p reorder-qfs
 spu_2p_SOURCES   = spu_2p.cpp
 endif

examples/reorder-qfs.c

@@ -0,0 +1,161 @@
+/* Copyright 2011 (c) Jostein R. Natvig <jostein.natvig@gmail.com> */
+
+
+
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <opm/core/grid.h>
+#include <opm/core/pressure/tpfa/ifs_tpfa.h>
+#include <opm/core/linalg/sparse_sys.h>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+
+#include <opm/core/linalg/call_umfpack.h>
+#include <opm/core/utility/cart_grid.h>
+
+#include <opm/core/transport/reorder/twophasetransport.h>
+
+
+struct Rock {
+    double *perm;
+    double *poro;
+};
+
+static void
+destroy_rock(struct Rock *rock)
+{
+    if (rock!=NULL)
+    {
+        free(rock->perm);
+        free(rock->poro);
+    }
+    free(rock);
+}
+
+static struct Rock*
+init_rock(int nc, int d)
+{
+    int i,j;
+    struct Rock *rock = malloc(sizeof *rock);
+
+    if (rock!=NULL)
+    {
+        rock->perm = malloc(nc*d*d*sizeof *rock->perm);
+        rock->poro = malloc(nc * sizeof *rock->poro);
+
+        if ((rock->perm==NULL) || (rock->poro==NULL))
+        {
+            destroy_rock(rock);
+            rock = NULL;
+        }
+        else
+        {
+            vector_zero(nc*d*d, rock->perm);
+            for (i=0; i<nc; ++i)
+            {
+                rock->poro[i] = 1.0;
+                for (j=0; j<d; ++j)
+                {
+                    rock->perm[j*(d+1)+i*d*d] = 1.0;
+                }
+            }
+        }
+    }
+
+    return rock; 
+}
+
+
+static void
+qfs(grid_t *g, struct Rock *rock, double *src, 
+    double *pressure, double *faceflux)
+{
+    size_t  alloc_sz, totconn;
+    double *htrans, *trans, *gpress;
+    struct ifs_tpfa_data *h;
+
+    fprintf(stderr, "QFS: Allocate space\n");
+    totconn = g->cell_facepos[ g->number_of_cells ];
+
+    alloc_sz  = totconn;            /* htrans */
+    alloc_sz += g->number_of_faces; /* trans */
+//    alloc_sz += g->number_of_cells; /* src */
+    alloc_sz += totconn;            /* gpress */
+
+    htrans = malloc(alloc_sz * sizeof *htrans);
+    trans  = htrans + totconn;
+//    src    = trans  + g->number_of_faces;
+    gpress = trans    + g->number_of_cells;
+
+    h = ifs_tpfa_construct(g);
+
+    fprintf(stderr, "QFS: Compute transmissibilities\n");
+    tpfa_htrans_compute(g, rock->perm, htrans);
+    tpfa_trans_compute (g, htrans    , trans );
+
+    vector_zero(g->number_of_cells, src   );
+    vector_zero(totconn           , gpress);
+
+    assert (g->number_of_cells > 1);
+
+    fprintf(stderr, "QFS: Assemble pressure system\n");
+    ifs_tpfa_assemble(g, trans, src, gpress, h);
+
+    fprintf(stderr, "QFS: Solve pressure system\n");
+    call_UMFPACK(h->A, h->b, h->x);
+
+    fprintf(stderr, "QFS: Recover pressure and flux\n");
+    ifs_tpfa_press_flux(g, trans, h, pressure, faceflux);
+
+    fprintf(stderr, "QFS: Release memory resources\n");
+    ifs_tpfa_destroy(h);
+    free(htrans);
+
+    fprintf(stderr, "QFS: Done\n");
+}
+
+void *compute_porevolume(struct UnstructuredGrid *g, double *poro, 
+    double *pv)
+{
+    int i;
+    for (i=0; i<g->number_of_cells; ++i)
+    {
+        pv[i] = g->cell_volumes[i]*poro[i];
+    }
+}
+
+int main(void)
+{
+    struct UnstructuredGrid *g = create_cart_grid_2d(2,1);
+    struct Rock *rock = init_rock(g->number_of_cells, g->dimensions);
+    double *sat   = malloc(g->number_of_cells *sizeof *sat);
+    double *press = malloc(g->number_of_cells *sizeof *press);
+    double *flux  = malloc(g->number_of_faces *sizeof *flux);
+    double *pv    = malloc(g->number_of_faces *sizeof *pv);
+    double *src    = malloc(g->number_of_faces *sizeof *src);
+
+    vector_zero(g->number_of_cells, sat);
+    vector_zero(g->number_of_cells, src);
+    src[0]=1;src[g->number_of_cells-1]=-1;
+
+    qfs(g, rock, src, press, flux);
+    
+    vector_write_stream(g->number_of_cells, press, stderr);
+
+    compute_porevolume(g, rock->poro, pv);
+
+    twophasetransport(pv, src, 1.0, g, flux, sat);
+
+    vector_write_stream(g->number_of_cells, sat, stderr);
+
+    free(pv);
+    free(src);
+    free(sat);
+    free(press);
+    free(flux);
+    destroy_rock(rock);
+    return 0;
+}