Merge OPMPressure into OPM core library. Unmodified.

src/cfs_tpfa.h

@@ -0,0 +1,128 @@
+/*
+  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_CFS_TPFA_HEADER_INCLUDED
+#define OPM_CFS_TPFA_HEADER_INCLUDED
+
+#include "grid.h"
+#include "flow_bc.h"
+#include "well.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct cfs_tpfa_impl;
+struct CSRMatrix;
+struct compr_quantities;
+
+struct cfs_tpfa_data {
+    struct CSRMatrix     *A;
+    double               *b;
+    double               *x;
+
+    struct cfs_tpfa_impl *pimpl;
+};
+
+
+struct cfs_tpfa_data *
+cfs_tpfa_construct(grid_t *G, well_t *W, int nphases);
+
+void
+cfs_tpfa_assemble(grid_t                  *G,
+                  double                   dt,
+                  well_t                  *W,
+                  flowbc_t                *bc,
+                  const double            *src,
+                  struct compr_quantities *cq,
+                  const double            *trans,
+                  const double            *gravcap_f,
+                  well_control_t          *wctrl,
+                  struct completion_data  *wdata,
+                  const double            *cpress0,
+                  const double            *porevol,
+                  struct cfs_tpfa_data    *h);
+
+void
+cfs_tpfa_press_flux(grid_t                 *G,
+                    flowbc_t               *bc,
+                    well_t                 *W,
+                    int                     np,
+                    const double           *trans,
+                    const double           *pmobf,
+                    const double           *gravcap_f,
+                    struct completion_data *wdata,
+                    struct cfs_tpfa_data   *h,
+                    double                 *cpress,
+                    double                 *fflux,
+                    double                 *wpress,
+                    double                 *wflux);
+
+void
+cfs_tpfa_fpress(grid_t               *G,
+                flowbc_t             *bc,
+                int                   np,
+                const double         *htrans,
+                const double         *pmobf,
+                const double         *gravcap_f,
+                struct cfs_tpfa_data *h,
+                const double         *cpress,
+                const double         *fflux,
+                double               *fpress);
+
+void
+cfs_tpfa_retrieve_masstrans(grid_t               *G,
+                            int                   np,
+                            struct cfs_tpfa_data *h,
+                            double               *masstrans_f);
+
+void
+cfs_tpfa_retrieve_gravtrans(grid_t               *G,
+                            int                   np,
+                            struct cfs_tpfa_data *h,
+                            double               *gravtrans_f);
+
+double
+cfs_tpfa_impes_maxtime(grid_t                  *G,
+                       struct compr_quantities *cq,
+                       const double            *trans,
+                       const double            *porevol,
+                       struct cfs_tpfa_data    *h,
+                       const double            *dpmobf,
+                       const double            *surf_dens,
+                       const double            *gravity);
+
+void
+cfs_tpfa_expl_mass_transport(grid_t                 *G,
+                             well_t                 *W,
+                             struct completion_data *wdata,
+                             int                     np,
+                             double                  dt,
+                             const double           *porevol,
+                             struct cfs_tpfa_data   *h,
+                             double                 *surf_vol);
+
+void
+cfs_tpfa_destroy(struct cfs_tpfa_data *h);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_CFS_TPFA_HEADER_INCLUDED */

src/cfs_tpfa_residual.h

@@ -0,0 +1,139 @@
+/*
+  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_CFS_TPFA_HEADER_INCLUDED
+#define OPM_CFS_TPFA_HEADER_INCLUDED
+
+#include "grid.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct cfs_tpfa_res_impl;
+struct CSRMatrix;
+struct compr_quantities_gen;
+struct compr_src;
+struct WellCompletions;
+struct WellControls;
+struct completion_data;
+
+struct cfs_tpfa_res_wells {
+    struct WellCompletions *conn;
+    struct WellControls    *ctrl;
+    struct completion_data *data;
+};
+
+struct cfs_tpfa_res_forces {
+    struct cfs_tpfa_res_wells *W  ;
+    struct compr_src          *src;
+};
+
+struct cfs_tpfa_res_data {
+    struct CSRMatrix         *J;
+    double                   *F;
+
+    struct cfs_tpfa_res_impl *pimpl;
+};
+
+
+struct cfs_tpfa_res_data *
+cfs_tpfa_res_construct(grid_t                 *G      ,
+                       struct WellCompletions *wconn  ,
+                       int                     nphases);
+
+void
+cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
+
+void
+cfs_tpfa_res_assemble(grid_t                      *G,
+                      double                       dt,
+                      struct cfs_tpfa_res_forces  *forces,
+                      const double                *zc,
+                      struct compr_quantities_gen *cq,
+                      const double                *trans,
+                      const double                *gravcap_f,
+                      const double                *cpress,
+                      const double                *wpress,
+                      const double                *porevol,
+                      struct cfs_tpfa_res_data    *h);
+
+void
+cfs_tpfa_res_flux(grid_t                     *G        ,
+                  struct cfs_tpfa_res_forces *forces   ,
+                  int                         np       ,
+                  const double               *trans    ,
+                  const double               *pmobf    ,
+                  const double               *gravcap_f,
+                  const double               *cpress   ,
+                  const double               *wpress   ,
+                  double                     *fflux    ,
+                  double                     *wflux    );
+
+void
+cfs_tpfa_res_fpress(grid_t                   *G,
+                    int                       np,
+                    const double             *htrans,
+                    const double             *pmobf,
+                    const double             *gravcap_f,
+                    struct cfs_tpfa_res_data *h,
+                    const double             *cpress,
+                    const double             *fflux,
+                    double                   *fpress);
+
+#if 0
+void
+cfs_tpfa_retrieve_masstrans(grid_t               *G,
+                            int                   np,
+                            struct cfs_tpfa_data *h,
+                            double               *masstrans_f);
+
+void
+cfs_tpfa_retrieve_gravtrans(grid_t               *G,
+                            int                   np,
+                            struct cfs_tpfa_data *h,
+                            double               *gravtrans_f);
+
+double
+cfs_tpfa_impes_maxtime(grid_t                  *G,
+                       struct compr_quantities *cq,
+                       const double            *trans,
+                       const double            *porevol,
+                       struct cfs_tpfa_data    *h,
+                       const double            *dpmobf,
+                       const double            *surf_dens,
+                       const double            *gravity);
+
+void
+cfs_tpfa_expl_mass_transport(grid_t                 *G,
+                             well_t                 *W,
+                             struct completion_data *wdata,
+                             int                     np,
+                             double                  dt,
+                             const double           *porevol,
+                             struct cfs_tpfa_data   *h,
+                             double                 *surf_vol);
+
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_CFS_TPFA_HEADER_INCLUDED */

src/compr_quant.h

@@ -0,0 +1,64 @@
+/*
+  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_COMPR_QUANT_HEADER_INCLUDED
+#define OPM_COMPR_QUANT_HEADER_INCLUDED
+
+#include <stddef.h>
+
+#include "grid.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct compr_quantities {
+    int           nphases;      /* Number of phases/components */
+
+    double *totcompr;     /* Total compressibility per cell */
+    double *voldiscr;     /* Volume discrepancy per cell */
+    double *Ac;           /* RB^{-1} per cell */
+    double *Af;           /* RB^{-1} per face */
+    double *phasemobf;    /* Phase mobility per face */
+};
+
+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);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_COMPR_QUANT_HEADER_INCLUDED */

src/compr_quant_general.h

@@ -0,0 +1,49 @@
+/*
+  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_COMPR_QUANT_HEADER_INCLUDED
+#define OPM_COMPR_QUANT_HEADER_INCLUDED
+
+#include <stddef.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct compr_quantities_gen {
+    int     nphases;      /* Number of phases/components */
+
+    double *Ac;           /* RB^{-1} per cell */
+    double *dAc;          /* d/dp (RB^{-1}) per cell */
+    double *Af;           /* RB^{-1} per face */
+    double *phasemobf;    /* Phase mobility per face */
+    double *voldiscr;     /* Volume discrepancy per cell */
+};
+
+struct compr_quantities_gen *
+compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
+
+void
+compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_COMPR_QUANT_HEADER_INCLUDED */

src/compr_source.h

@@ -0,0 +1,76 @@
+/*===========================================================================
+//
+// File: compr_source.h
+//
+// Created: 2011-10-19 19:14:30+0200
+//
+// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
+//          Jostein R. Natvig     <Jostein.R.Natvig@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  Copyright 2011 SINTEF ICT, Applied Mathematics.
+  Copyright 2011 Statoil ASA.
+
+  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_COMPR_SOURCE_H_HEADER
+#define OPM_COMPR_SOURCE_H_HEADER
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct compr_src {
+    int     nsrc;
+    int     cpty;
+
+    int     nphases;
+
+    int    *cell;
+    double *flux;
+    double *saturation;
+};
+
+
+struct compr_src *
+compr_src_allocate(int np, int nsrc);
+
+void
+compr_src_deallocate(struct compr_src *src);
+
+int
+append_compr_source_term(int               c  ,
+                         int               np ,
+                         double            v  ,
+                         const double     *sat,
+                         struct compr_src *src);
+
+void
+clear_compr_source_term(struct compr_src *src);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_COMPR_SOURCE_H_HEADER */

src/dfs.c

@@ -0,0 +1,118 @@
+/*
+  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/>.
+*/
+
+
+#include <assert.h>
+
+#include "dfs.h"
+
+/* 
+ * Assign color (nonnegative number) to each connected component of graph 
+ */
+void dfs (int size, int *ia, int *ja, int *ncolors, int *color, int* work)
+{
+    int i, c;
+    enum {UNVISITED = -1, VISITED = -2};
+    int *stack  = work;
+    int *count  = work + size;
+    int *bottom = stack;
+
+    *ncolors = 0; /* colors are nonnegative */
+  
+    for (i=0; i<size; ++i) {
+        color [i] = UNVISITED;
+        count [i] = ia[i+1]-ia[i];
+    }
+
+    /* Push seeds on stack */
+    for (i=0; i<size; ++i) {
+        if(color[i] >= 0) {  /* FINISHED */
+            continue;
+        }
+
+        *stack++ = i; /* push i */
+        color[i] = VISITED;
+
+        while ( stack != bottom ) {
+            c = *(stack-1); /* peek */
+            
+            if (count[c] > 0){
+                int child = ja[ia[c] + count[c]-1];
+                count[c]--;
+                
+                if (color[child] == UNVISITED) {
+                    *stack++ = child;
+                    color[c] = VISITED;
+                }
+   
+            } else {
+                color[c] = *ncolors;
+                --stack; /* pop c */
+            }
+        }
+        ++*ncolors;
+    }
+}
+
+
+
+
+#if defined(TEST) && TEST
+#include <stdlib.h>
+#include <stdio.h>
+
+
+/* Test code.  Reads a sparse matrix from a file specified on the   */
+/* command line with file format "m n\n i1 j1 v1\n i2 j2 v2\n ...". */
+/* Computes reorder sequence                                        */
+int main (int argc, char *argv [])
+{
+    int      *color, *work;
+    int       j, ncolors;
+#if 0
+    int size = 8;
+    int ia[] = {0, 1, 2, 4, 5, 5, 7, 7, 8};
+    int ja[] = {1, 2, 0, 3, 4, 5, 6, 6};
+#else
+    int size = 3;
+    int ia[] = {0,2,5,7};
+    int ja[] = {0,1,1,0,2,2,1};
+#endif
+
+
+    color = malloc (size * sizeof *color);
+    work  = malloc (2*size * sizeof *work);
+    dfs(size, ia, ja, &ncolors, color, work);
+
+
+    fprintf(stderr, "ncolors = %d\n", ncolors);
+    for (j=0; j<size; ++j) {
+        fprintf(stderr, "%d\n", color[j]);   
+    }
+ 
+
+    free (color);
+    free (work);
+    return 0;
+}
+#endif
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

src/dfs.h

@@ -0,0 +1,35 @@
+/*
+  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_DFS_HEADER_INCLUDED
+#define OPM_DFS_HEADER_INCLUDED
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+void dfs (int size, int *ia, int *ja, int *ncolors, int *color, int* work);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

src/flow_bc.c

@@ -0,0 +1,70 @@
+/*
+  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/>.
+*/
+
+#include <stdlib.h>
+
+#include "flow_bc.h"
+
+
+/* Create structure to hold flow boundary conditions for 'nf' faces.
+ *
+ * Return fully allocated structure or NULL in case of allocation
+ * failure. */
+/* ---------------------------------------------------------------------- */
+flowbc_t *
+allocate_flowbc(size_t nf)
+/* ---------------------------------------------------------------------- */
+{
+    size_t    i;
+    flowbc_t *new;
+
+    new = malloc(1 * sizeof *new);
+    if (new != NULL) {
+        new->type  = malloc(nf * sizeof *new->type);
+        new->bcval = malloc(nf * sizeof *new->bcval);
+
+        if ((new->type == NULL) || (new->bcval == NULL)) {
+            deallocate_flowbc(new);
+            new = NULL;
+        } else {
+            for (i = 0; i < nf; i++) {
+                new->type [i] = UNSET;
+                new->bcval[i] = 0.0;
+            }
+        }
+    }
+
+    return new;
+}
+
+
+/* Release memory resources for dynamically allocated flow boundary
+ * condition structure. */
+/* ---------------------------------------------------------------------- */
+void
+deallocate_flowbc(flowbc_t *fbc)
+/* ---------------------------------------------------------------------- */
+{
+    if (fbc != NULL) {
+        free(fbc->bcval);
+        free(fbc->type );
+    }
+
+    free(fbc);
+}

src/flow_bc.h

@@ -0,0 +1,49 @@
+/*
+  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_FLOW_BC_HEADER_INCLUDED
+#define OPM_FLOW_BC_HEADER_INCLUDED
+
+#include <stddef.h>
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+enum flowbc_type { UNSET, PRESSURE, FLUX };
+
+typedef struct {
+    enum flowbc_type *type;
+    double           *bcval;
+} flowbc_t;
+
+
+flowbc_t *
+allocate_flowbc(size_t nf);
+
+void
+deallocate_flowbc(flowbc_t *fbc);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_FLOW_BC_HEADER_INCLUDED */

src/ifs_tpfa.c

@@ -0,0 +1,273 @@
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#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;
+}
+
+
+/* ---------------------------------------------------------------------- */
+/* fgrav = accumarray(cf(j), grav(j).*sgn(j), [nf, 1]) */
+/* ---------------------------------------------------------------------- */
+static void
+compute_grav_term(grid_t *G, const double *gpress,
+                  double *fgrav)
+/* ---------------------------------------------------------------------- */
+{
+    int    c, i, f, c1, c2;
+    double s;
+
+    vector_zero(G->number_of_faces, fgrav);
+
+    for (c = i = 0; c < G->number_of_cells; c++) {
+        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;
+
+            if ((c1 >= 0) && (c2 >= 0)) {
+                fgrav[f] += s * gpress[i];
+            }
+        }
+    }
+}
+
+
+
+/* ======================================================================
+ * 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;
+
+    double s;
+
+    csrmatrix_zero(         h->A);
+    vector_zero   (h->A->m, h->b);
+
+    compute_grav_term(G, gpress, 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];
+
+            s  = 2.0*(c1 == c) - 1.0;
+            c2 = (c1 == c) ? c2 : c1;
+
+            h->b[c] -= trans[f] * (s * h->pimpl->fgrav[f]);
+
+            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;
+    double dh;
+
+    /* 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)) {
+            dh = cpress[c1] - cpress[c2] + h->pimpl->fgrav[f];
+            fflux[f] = trans[f] * dh;
+        } 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);
+}

src/ifs_tpfa.h

@@ -0,0 +1,65 @@
+/*
+  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"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+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);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_IFS_TPFA_HEADER_INCLUDED */

src/mimetic.c

@@ -0,0 +1,252 @@
+/*
+  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/>.
+*/
+
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+
+#include "blas_lapack.h"
+#include "mimetic.h"
+
+/* ------------------------------------------------------------------ */
+void
+mim_ip_simple_all(int ncells, int d, int max_nconn,
+                  int *pconn, int *conn,
+                  int *fneighbour, double *fcentroid, double *fnormal,
+                  double *farea, double *ccentroid, double *cvol,
+                  double *perm, double *Binv)
+/* ------------------------------------------------------------------ */
+{
+    int i, j, c, f, nconn, fpos2, lwork;
+
+    double *C, *N, *A, *work, s;
+
+    double cc[3] = { 0.0 };     /* No more than 3 space dimensions */
+
+    lwork = 64 * (max_nconn * d);                 /* 64 from ILAENV() */
+    C     = malloc((max_nconn * d) * sizeof *C);
+    N     = malloc((max_nconn * d) * sizeof *N);
+    A     = malloc(max_nconn       * sizeof *A);
+    work  = malloc(lwork           * sizeof *work);
+
+    if ((C != NULL) && (N != NULL) && (A != NULL) && (work != NULL)) {
+        fpos2 = 0;
+
+        for (c = 0; c < ncells; c++) {
+            for (j = 0; j < d; j++) {
+                cc[j] = ccentroid[j + c*d];
+            }
+
+            nconn = pconn[c + 1] - pconn[c];
+
+            for (i = 0; i < nconn; i++) {
+                f = conn[pconn[c] + i];
+                s = 2.0*(fneighbour[2 * f] == c) - 1.0;
+
+                A[i] = farea[f];
+
+                for (j = 0; j < d; j++) {
+                    C[i + j*nconn] = fcentroid  [j + f*d] - cc[j];
+                    N[i + j*nconn] = s * fnormal[j + f*d];
+                }
+            }
+
+            mim_ip_simple(nconn, nconn, d, cvol[c], &perm[c * d * d],
+                          C, A, N, &Binv[fpos2], work, lwork);
+
+            fpos2 += nconn * nconn;
+        }
+    }
+
+    free(work);  free(A);  free(N);  free(C);
+}
+
+
+/* ------------------------------------------------------------------ */
+void
+mim_ip_simple(int nf, int nconn, int d,
+              double v, double *K, double *C,
+              double *A, double *N,
+              double *Binv,
+              double *work, int lwork)
+/* ------------------------------------------------------------------ */
+{
+    mim_ip_span_nullspace(nf, nconn, d,    C, A, Binv, work, lwork);
+    mim_ip_linpress_exact(nf, nconn, d, v, K, N, Binv, work, lwork);
+}
+
+
+/* ------------------------------------------------------------------ */
+void
+mim_ip_span_nullspace(int nf, int nconn, int d,
+                      double *C,
+                      double *A,
+                      double *X,
+                      double *work, int nwork)
+/* ------------------------------------------------------------------ */
+{
+    MAT_SIZE_T m, n, k, ldC, ldX, info, lwork;
+
+    int    i, j;
+    double a1, a2;
+
+    double tau[3] = { 0.0 };  /* No more than 3 spatial dimensions */
+
+    /* Step 1) X(1:nf, 1:nf) <- I_{nf} */
+    for (j = 0; j < nf; j++) {
+        for (i = 0; i < nf; i++) {
+            X[i + j*nconn] = 0.0;
+        }
+        X[j * (nconn + 1)] = 1.0;
+    }
+
+    /* Step 2) C <- orth(A * C) */
+    for (j = 0; j < d; j++) {
+        for (i = 0; i < nf; i++) {
+            C[i + j*nf] *= A[i];
+        }
+    }
+
+    m = nf;  n = d;  ldC = nf;  k = d;  lwork = nwork;
+    dgeqrf_(&m, &n,     C, &ldC, tau, work, &lwork, &info);
+    dorgqr_(&m, &n, &k, C, &ldC, tau, work, &lwork, &info);
+
+    /* Step 3) X <- A * (X - C*C') * A */
+    ldX = nconn;
+    a1 = -1.0;  a2 = 1.0;
+    dsyrk_("Upper Triangular", "No Transpose",
+           &m, &n, &a1, C, &ldC, &a2, X, &ldX);
+    for (j = 0; j < nf; j++) {
+        for (i = 0; i <= j; i++) {
+            X[i + j*nconn] *= A[i] * A[j];
+        }
+    }
+
+    /* Account for DSYRK only assigning upper triangular part. */
+    for (j = 0; j < nf; j++) {
+        for (i = j + 1; i < nf; i++) {
+            X[i + j*nconn] = X[j + i*nconn];
+        }
+    }
+}
+
+
+/* ------------------------------------------------------------------ */
+void
+mim_ip_linpress_exact(int nf, int nconn, int d,
+                      double vol, double *K,
+                      double *N,
+                      double *Binv,
+                      double *work, int lwork)
+/* ------------------------------------------------------------------ */
+{
+    MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
+    int        i;
+    double     a1, a2, t;
+
+    assert (lwork >= d * nf);
+#if defined(NDEBUG)
+    /* Suppress warning about unused parameter. */
+    (void) lwork;
+#endif
+
+    t = 0.0;
+    for (i = 0; i < d; i++) {
+        t += K[i + i*d];
+    }
+
+    /* Step 4) T <- N*K */
+    m   = nf ;  n   = d  ;  k = d;
+    ld1 = nf ;  ld2 = d  ;
+    a1  = 1.0;  a2  = 0.0;
+    dgemm_("No Transpose", "No Transpose", &m, &n, &k,
+           &a1, N, &ld1, K, &ld2, &a2, work, &ld1);
+
+    /*  Step 5) Binv <- (N*K*N' + t*X) / vol */
+    a1 = 1.0     /      vol ;
+    a2 = 6.0 * t / (d * vol);
+    ldBinv = nconn;
+    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)_0 */
+/* ---------------------------------------------------------------------- */
+void
+mim_ip_mobility_update(int nc, const int *pconn, const double *totmob,
+                       const double *Binv0, 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] * Binv0[p2 + i];
+        }
+
+        p2 += n * n;
+    }
+}
+
+
+/* G <- \sum_i \rho_i f_i(s) * G_0 */
+/* ---------------------------------------------------------------------- */
+void
+mim_ip_density_update(int nc, const int *pconn, const double *omega,
+                      const double *gpress0, double *gpress)
+/* ---------------------------------------------------------------------- */
+{
+    int c, i;
+
+    for (c = i = 0; c < nc; c++) {
+        for (; i < pconn[c + 1]; i++) {
+            gpress[i] = omega[c] * gpress0[i];
+        }
+    }
+}

src/mimetic.h

@@ -0,0 +1,103 @@
+/*
+  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_MIMETIC_HEADER_INCLUDED
+#define OPM_MIMETIC_HEADER_INCLUDED
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void mim_ip_span_nullspace(int nf, int nconn, int d,
+                           double *C,
+                           double *A,
+                           double *X,
+                           double *work, int lwork);
+
+void mim_ip_linpress_exact(int nf, int nconn, int d,
+                           double vol, double *K,
+                           double *N,
+                           double *Binv,
+                           double *work, int lwork);
+
+void mim_ip_simple(int nf, int nconn, int d,
+                   double v, double *K, double *C,
+                   double *A, double *N,
+                   double *Binv,
+                   double *work, int lwork);
+
+
+/** Compute the mimetic inner products given a grid and cellwise
+ *  permeability tensors.
+ *
+ * @param ncells Number of cells in grid.
+ * @param d Number of space dimensions.
+ * @param max_ncf Maximum number of faces per cell.
+ * @param pconn Start indices in conn for each cell, plus end
+ *              marker. The size of pconn is (ncells + 1), and for a
+ *              cell i, [conn[pconn[i]], conn[pconn[i+1]]) is a
+ *              half-open interval containing the indices of faces
+ *              adjacent to i.
+ * @param conn Cell to face mapping. Size shall be equal to the sum of
+ *             ncf. See pconn for explanation.
+ * @param fneighbour Face to cell mapping. Its size shall be equal to
+ *                   the number of faces times 2. For each face, the
+ *                   two entries are either a cell number or -1
+ *                   (signifying the outer boundary). The face normal
+ *                   points out of the first cell and into the second.
+ * @param fcentroid Face centroids. Size shall be equal to the number
+ *                  of faces times d.
+ * @param fnormal Face normale. Size shall be equal to the number
+ *                of faces times d.
+ * @param farea Face areas.
+ * @param ccentroid Cell centroids. Size shall be ncells*d.
+ * @param cvol Cell volumes.
+ * @param perm Permeability. Size shall be ncells*d*d, storing a
+ *             d-by-d positive definite tensor per cell.
+ * @param[out] Binv This is where the inner product will be
+ *                  stored. Its size shall be equal to \f$\sum_i
+ * n_i^2\f$.
+ */
+void mim_ip_simple_all(int ncells, int d, int max_ncf,
+                       int *pconn, int *conn,
+                       int *fneighbour, double *fcentroid, double *fnormal,
+                       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)_0 */
+void
+mim_ip_mobility_update(int nc, const int *pconn, const double *totmob,
+                       const double *Binv0, double *Binv);
+
+/* G <- \sum_i \rho_i f_i(s) * G_0 */
+void
+mim_ip_density_update(int nc, const int *pconn, const double *omega,
+                      const double *gpress0, double *gpress);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* OPM_MIMETIC_HEADER_INCLUDED */

src/partition.c

@@ -0,0 +1,581 @@
+/*
+  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/>.
+*/
+
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "dfs.h"
+#include "partition.h"
+
+
+#define MAX(a,b) (((a) > (b)) ? (a) : (b))
+
+
+/* [cidx{1:ndims}] = ind2sub(size, idx) */
+/* ---------------------------------------------------------------------- */
+static void
+partition_coord_idx(int ndims, int idx, const int *size, int *cidx)
+/* ---------------------------------------------------------------------- */
+{
+    int i;
+
+    for (i = 0; i < ndims; i++) {
+        cidx[i]  = idx % size[i];
+        idx     /=       size[i];
+    }
+
+    assert (idx == 0);
+}
+
+
+/* sub2ind(size, cidx{1:ndims}) */
+/* ---------------------------------------------------------------------- */
+static int
+partition_lin_idx(int ndims, const int *size, const int *cidx)
+/* ---------------------------------------------------------------------- */
+{
+    int i, idx;
+
+    idx = cidx[ndims - 1];
+    for (i = ndims - 2; i >= 0; i--) {
+        idx = cidx[i] + size[i]*idx;
+    }
+
+    return idx;
+}
+
+
+/* ---------------------------------------------------------------------- */
+/* Load-balanced linear distribution.
+ *
+ * See Eric F. Van de Velde, Concurrent Scientific Computing,
+ * 1994, Springer Verlag, p. 54 (Sect. 2.3) for details.  */
+static void
+partition_loadbal_lin_dist(int ndims, const int *size, const int *nbins,
+                           int *idx)
+/* ---------------------------------------------------------------------- */
+{
+    int i, L, R, b1, b2;
+
+    for (i = 0; i < ndims; i++) {
+        L = size[i] / nbins[i]; /* # entities per bin */
+        R = size[i] % nbins[i]; /* # bins containing one extra entity */
+
+        b1 =  idx[i]      / (L + 1);
+        b2 = (idx[i] - R) /  L     ;
+
+        idx[i] = MAX(b1, b2);
+    }
+}
+
+
+/* Partition 'nc' fine-scale Cartesian indices 'idx' from a box of
+ * dimension 'fine_d' into a coarse-scale box of dimension 'coarse_d'.
+ *
+ * Store partition in vector 'p' (assumed to hold at least 'nc'
+ * slots).
+ *
+ * Allocates a tiny work array to hold 'ndims' ints.  Returns 'nc' if
+ * successful and -1 if unable to allocate the work array. */
+/* ---------------------------------------------------------------------- */
+int
+partition_unif_idx(int ndims, int nc,
+                   const int *fine_d, const int *coarse_d, const int *idx,
+                   int *p)
+/* ---------------------------------------------------------------------- */
+{
+    int c, ret, *ix;
+
+    ix = malloc(ndims * sizeof *ix);
+
+    if (ix != NULL) {
+        for (c = 0; c < nc; c++) {
+            partition_coord_idx(ndims, idx[c], fine_d, ix);
+
+            partition_loadbal_lin_dist(ndims, fine_d, coarse_d, ix);
+
+            p[c] = partition_lin_idx(ndims, coarse_d, ix);
+        }
+
+        ret = nc;
+    } else {
+        ret = -1;
+    }
+
+    free(ix);
+
+    return ret;
+}
+
+
+/* Renumber blocks to create contiguous block numbers from 0..n-1
+ * (in other words: remove empty coarse blocks).
+ *
+ * Returns maximum new block number if successful and -1 if not. */
+/* ---------------------------------------------------------------------- */
+int
+partition_compress(int n, int *p)
+/* ---------------------------------------------------------------------- */
+{
+    int ret, i, max, *compr;
+
+    max = -1;
+    for (i = 0; i < n; i++) {
+        assert (0 <= p[i]);     /* Only non-neg partitions (for now?). */
+        max = MAX(max, p[i]);
+    }
+
+    compr = malloc((max + 1) * sizeof *compr);
+
+    if (compr != NULL) {
+        for (i = 0; i < max + 1; i++) { compr[i] = 0; }
+        for (i = 0; i < n; i++)       { compr[p[i]]++; }
+
+        compr[0] = -1 + (compr[0] > 0);
+        for (i = 1; i <= max; i++) {
+            compr[i] = compr[i - 1] + (compr[i] > 0);
+        }
+
+        for (i = 0; i < n; i++) { p[i] = compr[p[i]]; }
+
+        ret = compr[max];
+    } else {
+        ret = -1;
+    }
+
+    free(compr);
+
+    return ret;
+}
+
+
+/* Free memory resources for block->cell map. */
+/* ---------------------------------------------------------------------- */
+void
+partition_deallocate_inverse(int *pi, int *inverse)
+/* ---------------------------------------------------------------------- */
+{
+    free(inverse);
+    free(pi);
+}
+
+
+/* Allocate memory for block->cell map (CSR representation).  Highest
+ * block number is 'max_bin'.  Grid contains 'nc' cells.
+ *
+ * Returns 'nc' (and sets CSR pointer pair (*pi, *inverse)) if
+ * successful, -1 and pointers to NULL if not. */
+/* ---------------------------------------------------------------------- */
+int
+partition_allocate_inverse(int nc, int max_bin,
+                           int **pi, int **inverse)
+/* ---------------------------------------------------------------------- */
+{
+    int nbin, ret, *ptr, *i;
+
+    nbin = max_bin + 1;
+
+    ptr  = malloc((nbin + 1) * sizeof *ptr);
+    i    = malloc(nc         * sizeof *i  );
+
+    if ((ptr == NULL) || (i == NULL)) {
+        partition_deallocate_inverse(ptr, i);
+
+        *pi      = NULL;
+        *inverse = NULL;
+
+        ret = 0;
+    } else {
+        *pi      = ptr;
+        *inverse = i;
+
+        ret = nc;
+    }
+
+    return ret;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static int
+max_block(int nc, const int *p)
+/* ---------------------------------------------------------------------- */
+{
+    int m, i;
+
+    m = -1;
+
+    for (i = 0; i < nc; i++) {
+        m = MAX(m, p[i]);
+    }
+
+    return m;
+}
+
+
+/* Invert cell->block mapping 'p' (partition vector) to create
+ * block->cell mapping (CSR representation, pointer pair (pi,inverse)). */
+/* ---------------------------------------------------------------------- */
+void
+partition_invert(int nc, const int *p, int *pi, int *inverse)
+/* ---------------------------------------------------------------------- */
+{
+    int nbin, b, i;
+
+    nbin = max_block(nc, p) + 1; /* Adjust for bin 0 */
+
+    /* Zero start pointers */
+    for (i = 0; i < nbin + 1; i++) {
+        pi[i] = 0;
+    }
+
+    /* Count elements per bin */
+    for (i = 0; i < nc; i++) { pi[ p[i] + 1 ]++; }
+
+    for (b = 1; b <= nbin; b++) {
+        pi[0] += pi[b];
+        pi[b]  = pi[0] - pi[b];
+    }
+
+    /* Insert bin elements whilst deriving start pointers */
+    for (i = 0; i < nc; i++) {
+        inverse[ pi[ p[i] + 1 ] ++ ] = i;
+    }
+
+    /* Assert basic sanity */
+    assert (pi[nbin] == nc);
+    pi[0] = 0;
+}
+
+
+/* Create local cell numbering, within the cell's block, for each
+ * global cell. */
+/* ---------------------------------------------------------------------- */
+void
+partition_localidx(int nbin, const int *pi, const int *inverse,
+                   int *localidx)
+/* ---------------------------------------------------------------------- */
+{
+    int b, i;
+
+    for (b = 0; b < nbin; b++) {
+        for (i = pi[b]; i < pi[b + 1]; i++) {
+            localidx[ inverse[i] ] = i - pi[b];
+        }
+    }
+}
+
+
+/* Release memory resources for internal cell-to-cell connectivity
+ * (CSR representation). */
+/* ---------------------------------------------------------------------- */
+static void
+partition_destroy_c2c(int *pc2c, int *c2c)
+/* ---------------------------------------------------------------------- */
+{
+    free(c2c);  free(pc2c);
+}
+
+
+/* Create symmetric cell-to-cell (internal) connectivity for domain
+ * containing 'nc' cells.  CSR representation (*pc2c,*c2c).
+ *
+ * Neighbourship 'neigh' is 2*nneigh array such that cell neigh[2*i+0]
+ * is connected to cell neigh[2*i+1] for all i=0:nneigh-1.
+ *
+ * Negative 'neigh' entries represent invalid cells (outside domain).
+ *
+ * Returns 'nc' (and sets pointer pair) if successful, 0 (and pointer
+ * pair to NULL) if not. */
+/* ---------------------------------------------------------------------- */
+static int
+partition_create_c2c(int nc, int nneigh, const int *neigh,
+                     int **pc2c, int **c2c)
+/* ---------------------------------------------------------------------- */
+{
+    int i, ret, c1, c2;
+
+    *pc2c = malloc((nc + 1) * sizeof **pc2c);
+
+    if (*pc2c != NULL) {
+        for (i = 0; i < nc + 1; i++) {
+            (*pc2c)[i] = 0;
+        }
+
+        for (i = 0; i < nneigh; i++) {
+            c1 = neigh[2*i + 0];
+            c2 = neigh[2*i + 1];
+
+            if ((c1 >= 0) && (c2 >= 0)) {
+                /* Symmetric Laplace matrix (undirected graph) */
+                (*pc2c)[ c1 + 1 ] ++;
+                (*pc2c)[ c2 + 1 ] ++;
+            }
+        }
+
+        for (i = 1; i <= nc; i++) {
+            (*pc2c)[i] += 1;    /* Self connection */
+
+            (*pc2c)[0] += (*pc2c)[i];
+            (*pc2c)[i]  = (*pc2c)[0] - (*pc2c)[i];
+        }
+
+        *c2c = malloc((*pc2c)[0] * sizeof **c2c);
+
+        if (*c2c != NULL) {
+            /* Self connections */
+            for (i = 0; i < nc; i++) {
+                (*c2c)[ (*pc2c)[i + 1] ++ ] = i;
+            }
+
+            for (i = 0; i < nneigh; i++) {
+                c1 = neigh[2*i + 0];
+                c2 = neigh[2*i + 1];
+
+                if ((c1 >= 0) && (c2 >= 0)) {
+                    /* Symmetric Laplace matrix (undirected graph) */
+                    (*c2c)[ (*pc2c)[ c1 + 1 ] ++ ] = c2;
+                    (*c2c)[ (*pc2c)[ c2 + 1 ] ++ ] = c1;
+                }
+            }
+
+            ret = nc;
+        } else {
+            free(*pc2c);
+            *pc2c = NULL;
+
+            ret = 0;
+        }
+    } else {
+        *c2c = NULL;
+
+        ret = 0;
+    }
+
+    return ret;
+}
+
+
+/* Release dfs() memory resources. */
+/* ---------------------------------------------------------------------- */
+static void
+deallocate_dfs_arrays(int *ia, int *ja, int *colour, int *work)
+/* ---------------------------------------------------------------------- */
+{
+    free(work);  free(colour);  free(ja);  free(ia);
+}
+
+
+/* Allocate dfs() memory resources to support graph containing 'n'
+ * nodes and (at most) 'nnz' total connections.  Return 'n' if
+ * successful (and set pointers) and 0 (and set pointers to NULL) if
+ * not. */
+/* ---------------------------------------------------------------------- */
+static int
+allocate_dfs_arrays(int n, int nnz,
+                    int **ia, int **ja, int **colour, int **work)
+/* ---------------------------------------------------------------------- */
+{
+    int ret;
+
+    *ia     = malloc((n + 1) * sizeof **ia    );
+    *ja     = malloc(nnz     * sizeof **ja    );
+    *colour = malloc(n       * sizeof **colour);
+    *work   = malloc(2 * n   * sizeof **work  );
+
+    if ((*ia     == NULL) || (*ja   == NULL) ||
+        (*colour == NULL) || (*work == NULL)) {
+        deallocate_dfs_arrays(*ia, *ja, *colour, *work);
+
+        *ia     = NULL;
+        *ja     = NULL;
+        *colour = NULL;
+        *work   = NULL;
+
+        ret = 0;
+    } else {
+        ret = n;
+    }
+
+    return ret;
+}
+
+
+/* Compute maximum number of cells (*max_blk_cells) and cell-to-cell
+ * connections (*max_blk_conn) over all blocks. */
+/* ---------------------------------------------------------------------- */
+static void
+count_block_conns(int nblk,
+                  const int *pb2c, const int *b2c, const int *pc2c,
+                  int *max_blk_cells, int *max_blk_conn)
+/* ---------------------------------------------------------------------- */
+{
+    int b, i, n_blk_conn;
+
+    *max_blk_cells = 0;
+    *max_blk_conn  = 0;
+
+    i = 0;                      /* == pb2c[0] */
+    for (b = 0; b < nblk; b++) {
+        n_blk_conn = 0;
+
+        for (; i < pb2c[b + 1]; i++) {
+            n_blk_conn += pc2c[b2c[i] + 1] - pc2c[b2c[i]];
+        }
+
+        *max_blk_cells = MAX(*max_blk_cells, pb2c[b + 1] - pb2c[b]);
+        *max_blk_conn  = MAX(*max_blk_conn , n_blk_conn);
+    }
+}
+
+
+/* Create block-internal (symmetric) connectivity graph (CSR
+ * representation ia,ja) for connected component labelling (used in
+ * splitting disconnected blocks). */
+/* ---------------------------------------------------------------------- */
+static void
+create_block_conns(int        b   ,
+                   const int *p   , const int *loc,
+                   const int *pb2c, const int *b2c,
+                   const int *pc2c, const int *c2c,
+                   int       *ia  , int       *ja )
+/* ---------------------------------------------------------------------- */
+{
+    int nc, c, i, j;
+
+    nc = pb2c[b + 1] - pb2c[b];
+
+    /* Clear start pointers */
+    for (i = 0; i < nc + 1; i++) {
+        ia[i] = 0;
+    }
+
+    for (i = pb2c[b]; i < pb2c[b + 1]; i++) {
+        c = b2c[i];   assert (loc[c] == i - pb2c[b]);
+
+        /* Self connections inserted in partition_create_c2c()) */
+        for (j = pc2c[c]; j < pc2c[c + 1]; j++) {
+            if (p[c2c[j]] == b) {
+                /* Connection internal to block 'b'.  Add */
+                ia[loc[c] + 1] ++;
+            }
+        }
+    }
+
+    assert (ia[0] == 0);
+
+    for (i = 1; i <= nc; i++) {
+        ia[0] += ia[i];
+        ia[i]  = ia[0] - ia[i];
+    }
+
+    for (i = pb2c[b]; i < pb2c[b + 1]; i++) {
+        c = b2c[i];
+
+        /* Create connections (self conn automatic) */
+        for (j = pc2c[c]; j < pc2c[c + 1]; j++) {
+            if (p[c2c[j]] == b) {
+                ja[ ia[loc[c] + 1] ++ ] = loc[c2c[j]];
+            }
+        }
+    }
+
+    ia[0] = 0;
+}
+
+
+/* Split disconnected coarse blocks.  Preserve block numbering where
+ * possible.
+ *
+ * Neighbourship definition 'neigh' is pointer to 2*nneigh array such
+ * that cell neigh[2*i+0] is connected to cell neigh[2*i+1] for all
+ * i=0:nneigh-1.  Negative entries in 'neigh' represent invalid cells
+ * (outside domain).
+ *
+ * Returns number of new blocks (0 if all blocks internally connected)
+ * if successful and -1 otherwise. */
+/* ---------------------------------------------------------------------- */
+int
+partition_split_disconnected(int nc, int nneigh, const int *neigh, int *p)
+/* ---------------------------------------------------------------------- */
+{
+    int inv_ok, c2c_ok, dfs_ok;
+    int i, b, ret, maxblk, ncolour, max_blk_cells, max_blk_conn;
+
+    int *pb2c, *b2c, *loc, *pc2c, *c2c;
+    int *ia, *ja, *colour, *work;
+
+    maxblk = max_block(nc, p);
+
+    inv_ok = partition_allocate_inverse(nc, maxblk, &pb2c, &b2c);
+    c2c_ok = partition_create_c2c(nc, nneigh, neigh, &pc2c, &c2c);
+    loc    = malloc(nc * sizeof *loc);
+
+    if (inv_ok && c2c_ok && (loc != NULL)) {
+        partition_invert(nc, p, pb2c, b2c);
+        partition_localidx(maxblk + 1, pb2c, b2c, loc);
+
+        count_block_conns(maxblk + 1, pb2c, b2c, pc2c,
+                          &max_blk_cells, &max_blk_conn);
+
+        dfs_ok = allocate_dfs_arrays(max_blk_cells, max_blk_conn,
+                                     &ia, &ja, &colour, &work);
+        if (dfs_ok) {
+            /* Target acquired.  Fire. */
+            ret = 0;
+
+            for (b = 0; b < maxblk + 1; b++) {
+                create_block_conns(b, p, loc, pb2c, b2c, pc2c, c2c, ia, ja);
+
+                dfs(pb2c[b + 1] - pb2c[b], ia, ja, &ncolour, colour, work);
+
+                if (ncolour > 1) {
+                    /* Block contains more than one component.  Assign
+                     * new block numbers for cells in components
+                     * 1:ncomp-1. */
+                    for (i = pb2c[b]; i < pb2c[b + 1]; i++) {
+                        if (colour[i - pb2c[b]] > 0) {
+                            p[b2c[i]] = maxblk + ret + colour[i - pb2c[b]];
+                        }
+                    }
+
+                    ret += ncolour - 1;
+                }
+            }
+        } else {
+            ret = -1;
+        }
+
+        deallocate_dfs_arrays(ia, ja, colour, work);
+    } else {
+        ret = -1;
+    }
+
+    free(loc);
+    partition_destroy_c2c(pc2c, c2c);
+    partition_deallocate_inverse(pb2c, b2c);
+
+    return ret;
+}
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

src/partition.h

@@ -0,0 +1,62 @@
+/*
+  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_PARTITION_HEADER_INCLUDED
+#define OPM_PARTITION_HEADER_INCLUDED
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int
+partition_unif_idx(int ndims, int nc,
+                   const int *fine_d,
+                   const int *coarse_d,
+                   const int *idx,
+                   int *p);
+
+int
+partition_compress(int n, int *p);
+
+
+int
+partition_allocate_inverse(int nc, int max_blk,
+                           int **pi, int **inverse);
+
+void
+partition_deallocate_inverse(int *pi, int *inverse);
+
+void
+partition_invert(int nc, const int *p,
+                 int *pi, int *inverse);
+
+void
+partition_localidx(int nblk, const int *pi, const int *inverse,
+                   int *localidx);
+
+
+int
+partition_split_disconnected(int nc, int nneigh, const int *neigh,
+                             int *p);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* PARTITION_H_INLCUDED */

src/sparse_sys.c

@@ -0,0 +1,263 @@
+/*
+  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/>.
+*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "sparse_sys.h"
+
+
+/* ---------------------------------------------------------------------- */
+struct CSRMatrix *
+csrmatrix_new_count_nnz(size_t m)
+/* ---------------------------------------------------------------------- */
+{
+    size_t            i;
+    struct CSRMatrix *new;
+
+    assert (m > 0);
+
+    new = malloc(1 * sizeof *new);
+    if (new != NULL) {
+        new->ia = malloc((m + 1) * sizeof *new->ia);
+
+        if (new->ia != NULL) {
+            for (i = 0; i < m + 1; i++) { new->ia[i] = 0; }
+
+            new->m   = m;
+            new->nnz = 0;
+
+            new->ja  = NULL;
+            new->sa  = NULL;
+        } else {
+            csrmatrix_delete(new);
+            new = NULL;
+        }
+    }
+
+    return new;
+}
+
+
+/* Allocate CSR matrix, known nnz.  Allocation only.  Caller must
+ * build sparsity structure before using in global assembly.
+ *
+ * Returns fully allocated structure if successful and NULL otherwise. */
+/* ---------------------------------------------------------------------- */
+struct CSRMatrix *
+csrmatrix_new_known_nnz(size_t m, size_t nnz)
+/* ---------------------------------------------------------------------- */
+{
+    struct CSRMatrix *new;
+
+    new = malloc(1 * sizeof *new);
+
+    if (new != NULL) {
+        new->ia = malloc((m + 1) * sizeof *new->ia);
+        new->ja = malloc(nnz     * sizeof *new->ja);
+        new->sa = malloc(nnz     * sizeof *new->sa);
+
+        if ((new->ia == NULL) || (new->ja == NULL) || (new->sa == NULL)) {
+            csrmatrix_delete(new);
+            new = NULL;
+        } else {
+            new->m   = m;
+            new->nnz = nnz;
+        }
+    }
+
+    return new;
+}
+
+
+/* ---------------------------------------------------------------------- */
+size_t
+csrmatrix_new_elms_pushback(struct CSRMatrix *A)
+/* ---------------------------------------------------------------------- */
+{
+    size_t i;
+
+    assert (A->ia[0] == 0);     /* Elems for row 'i' in bin i+1 ... */
+
+    for (i = 1; i <= A->m; i++) {
+        A->ia[0] += A->ia[i];
+        A->ia[i]  = A->ia[0] - A->ia[i];
+    }
+
+    A->nnz = A->ia[0];
+    assert (A->nnz > 0);        /* Else not a real system. */
+
+    A->ia[0] = 0;
+
+    A->ja = malloc(A->nnz * sizeof *A->ja);
+    A->sa = malloc(A->nnz * sizeof *A->sa);
+
+    if ((A->ja == NULL) || (A->sa == NULL)) {
+        free(A->sa);   A->sa = NULL;
+        free(A->ja);   A->ja = NULL;
+
+        A->nnz = 0;
+    }
+
+    return A->nnz;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static int
+cmp_row_elems(const void *a0, const void *b0)
+/* ---------------------------------------------------------------------- */
+{
+    return *(const int * const)a0 - *(const int * const)b0;
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+csrmatrix_sortrows(struct CSRMatrix *A)
+/* ---------------------------------------------------------------------- */
+{
+    size_t i;
+
+    /* O(A->nnz * log(average nnz per row)) \approx O(A->nnz) */
+    for (i = 0; i < A->m; i++) {
+        qsort(A->ja        + A->ia[i] ,
+              A->ia[i + 1] - A->ia[i] ,
+              sizeof A->ja  [A->ia[i]],
+              cmp_row_elems);
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+size_t
+csrmatrix_elm_index(int i, int j, const struct CSRMatrix *A)
+/* ---------------------------------------------------------------------- */
+{
+    int *p;
+
+    p = bsearch(&j, A->ja + A->ia[i], A->ia[i + 1] - A->ia[i],
+                sizeof A->ja[A->ia[i]], cmp_row_elems);
+
+    assert (p != NULL);
+
+    return p - A->ja;
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+csrmatrix_delete(struct CSRMatrix *A)
+/* ---------------------------------------------------------------------- */
+{
+    if (A != NULL) {
+        free(A->sa);
+        free(A->ja);
+        free(A->ia);
+    }
+
+    free(A);
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+csrmatrix_zero(struct CSRMatrix *A)
+/* ---------------------------------------------------------------------- */
+{
+    vector_zero(A->nnz, A->sa);
+}
+
+
+/* ---------------------------------------------------------------------- */
+/* 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; }
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+csrmatrix_write(const struct CSRMatrix *A, const char *fn)
+/* ---------------------------------------------------------------------- */
+{
+    FILE *fp;
+
+    fp = fopen(fn, "wt");
+
+    if (fp != NULL) {
+        csrmatrix_write_stream(A, fp);
+    }
+
+    fclose(fp);
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+csrmatrix_write_stream(const struct CSRMatrix *A, FILE *fp)
+/* ---------------------------------------------------------------------- */
+{
+    size_t  i, j;
+    for (i = j = 0; i < A->m; i++) {
+        for (; j < (size_t) (A->ia[i + 1]); j++) {
+            fprintf(fp, "%lu %lu %26.18e\n",
+                    (unsigned long) (i + 1),
+                    (unsigned long) (A->ja[j] + 1),
+                    A->sa[j]);
+        }
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+vector_write(size_t n, const double *v, const char *fn)
+/* ---------------------------------------------------------------------- */
+{
+    FILE *fp;
+
+    fp = fopen(fn, "wt");
+
+    if (fp != NULL) {
+        vector_write_stream(n, v, fp);
+    }
+
+    fclose(fp);
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+vector_write_stream(size_t n, const double *v, FILE *fp)
+/* ---------------------------------------------------------------------- */
+{
+    size_t i;
+
+    for (i = 0; i < n; i++) {
+        fprintf(fp, "%26.18e\n", v[i]);
+    }
+}

src/sparse_sys.h

@@ -0,0 +1,87 @@
+/*
+  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_SPARSE_SYS_HEADER_INCLUDED
+#define OPM_SPARSE_SYS_HEADER_INCLUDED
+
+#include <stddef.h>
+#include <stdio.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+struct CSRMatrix
+{
+    size_t      m;
+    size_t      nnz;
+
+    int        *ia;
+    int        *ja;
+
+    double     *sa;
+};
+
+
+
+struct CSRMatrix *
+csrmatrix_new_count_nnz(size_t m);
+
+struct CSRMatrix *
+csrmatrix_new_known_nnz(size_t m, size_t nnz);
+
+size_t
+csrmatrix_new_elms_pushback(struct CSRMatrix *A);
+
+size_t
+csrmatrix_elm_index(int i, int j, const struct CSRMatrix *A);
+
+void
+csrmatrix_sortrows(struct CSRMatrix *A);
+
+void
+csrmatrix_delete(struct CSRMatrix *A);
+
+void
+csrmatrix_zero(struct CSRMatrix *A);
+
+/* ---------------------------------------------------------------------- */
+/* v = zeros([n, 1]) */
+/* ---------------------------------------------------------------------- */
+void
+vector_zero(size_t n, double *v);
+
+void
+csrmatrix_write(const struct CSRMatrix *A, const char *fn);
+
+void
+csrmatrix_write_stream(const struct CSRMatrix *A, FILE *fp);
+
+void
+vector_write(size_t n, const double *v, const char *fn);
+
+void
+vector_write_stream(size_t n, const double *v, FILE *fp);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* OPM_SPARSE_SYS_HEADER_INCLUDED */

src/well.h

@@ -0,0 +1,67 @@
+/*
+  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_WELL_HEADER_INCLUDED
+#define OPM_WELL_HEADER_INCLUDED
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+enum well_type    { INJECTOR, PRODUCER };
+enum well_control { BHP     , RATE     };
+
+struct WellCompletions {
+    int  number_of_wells;
+    int *well_connpos;
+    int *well_cells;
+};
+
+struct WellControls {
+    enum well_type    *type;
+    enum well_control *ctrl;
+    double            *target;
+};
+
+struct completion_data {
+    double *WI;                 /* Productivity index */
+    double *gpot;               /* Gravity potential */
+    double *A;                  /* RB^{-1} */
+    double *phasemob;           /* Phase mobility */
+};
+
+typedef struct WellCompletions well_t;
+typedef struct WellControls    well_control_t;
+
+int
+allocate_cell_wells(int nc, well_t *W, int **cwpos, int **cwells);
+
+void
+deallocate_cell_wells(int *cvpos, int *cwells);
+
+void
+derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* OPM_WELL_HEADER_INCLUDED */