Add bare-bones implementation of reordering algorithm.

Lightly tested.

Makefile.am

@@ -73,8 +73,14 @@ opm/core/pressure/mimetic/mimetic.c                     \
 opm/core/pressure/mimetic/hybsys_global.c               \
 opm/core/pressure/mimetic/hybsys.c                      \
 opm/core/transport/spu_explicit.c                       \
-opm/core/transport/spu_implicit.c                       \
-opm/core/transport/transport_source.c
+opm/core/transport/spu_implicit.c		\
+opm/core/transport/transport_source.c		\
+opm/core/transport/reorder/reordersequence.c	\
+opm/core/transport/reorder/twophase.c		\
+opm/core/transport/reorder/nlsolvers.c		\
+opm/core/transport/reorder/tarjan.c		\
+opm/core/transport/reorder/twophasetransport.c
+
 
 nobase_include_HEADERS =                                \
 opm/core/eclipse/CornerpointChopper.hpp                 \
@@ -177,7 +183,12 @@ opm/core/transport/spu_explicit.h                       \
 opm/core/transport/JacobianSystem.hpp                   \
 opm/core/transport/spu_implicit.h                       \
 opm/core/transport/ImplicitAssembly.hpp                 \
-opm/core/transport/transport_source.h
+opm/core/transport/transport_source.h		\
+opm/core/transport/reorder/nlsolvers.h		\
+opm/core/transport/reorder/reordersequence.h	\
+opm/core/transport/reorder/twophase.h		\
+opm/core/transport/reorder/tarjan.h		\
+opm/core/transport/reorder/twophasetransport.h
 
 if UMFPACK
 libopmcore_la_SOURCES  += opm/core/linalg/call_umfpack.c

opm/core/transport/reorder/nlsolvers.c

@@ -0,0 +1,263 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+
+#include "nlsolvers.h"
+
+
+#ifdef MEXFILE
+#include <mex.h>
+extern int interrupted;
+#define print   mexPrintf
+#define malloc  mxMalloc
+#define calloc  mxCalloc
+#define realloc mxRealloc
+#define free    mxFree
+#else
+#define print printf
+#endif
+
+
+static const char no_root_str[]=
+  "  In %s:\n"
+  "  With G(0) =% 5f, G(1) =% 5f, G(s) cannot have a zero in [0,1]!\n";
+
+static const double EPS = 1e-14;
+
+
+/* Start with bracket [0,1] with G(0)*G(1)<0.  */
+/*---------------------------------------------------------------------------*/
+double
+ridder (double (*G)(double, void*), void *data, double tol, int maxit, int *iterations)
+/*---------------------------------------------------------------------------*/
+{
+  *iterations = 0;
+
+  double s2 = ((double *) data)[2];
+  double G2 = G(s2, data);
+  if (fabs(G2)<tol) return s2;
+
+  /* Initial guess is interval [0,1] of course */
+  double s0=0.0;
+  double G0=G(s0, data);
+  if (fabs(G0)<EPS) return s0;
+
+  double s1=1.0;
+  double G1=G(s1, data);
+  if (fabs(G1)<EPS) return s1;
+
+  if (G0*G1 > 0)
+  {
+    print(no_root_str, "ridder", G0, G1);
+    return -1.0;
+  }
+
+  if (G0>G1)
+  {
+    double swap;
+    swap = G0;
+    G0   = G1;
+    G1   = swap;
+  }
+
+  double s3=0;
+  double G3=10;
+
+  int it = 0;
+  while ( (fabs(G3) > tol)  &&  (it++ < maxit))
+  {
+    /* find zero crossing of line segment [(s0,G0), (s1,G1)] */
+    double root = sqrt(G2*G2-G0*G1);
+    if (fabs(root)<EPS)
+      return -1.0; /* Hmmm */
+
+    double sgn  = G0>G1 ? 1.0 : -1.0;
+    s3          = s2 + ( s2-s0 )*sgn*G2/root;
+    G3          = G(s3, data);
+
+
+    /* if     G2*G3<0  */
+    if (G2*G3 <= 0.0)
+    {
+      if (G2 > G3)
+      {
+	s0 = s3;
+	G0 = G3;
+	s1 = s2;
+	G1 = G2;
+      }
+      else
+      {
+	s0 = s2;
+	G0 = G2;
+	s1 = s3;
+	G1 = G3;
+
+      }
+
+    }
+    else if (G0*G3 <= 0.0)
+    {
+      s1 = s3;
+      G1 = G3;
+    }
+    else if (G1*G3 <= 0.0)
+    {
+      s0 = s3;
+      G0 = G3;
+    }
+    else
+    {
+      print(
+	    "In ridder:\n"
+	    "G0=%10.10f, G1=%10.10f, G3=%10.10f\n",
+	    G0, G1, G3
+	    );
+      getchar();
+    }
+    s2   = 0.5*(s0+s1);
+    G2   = G(s2, data);
+  }
+
+  *iterations = it;
+  return s3;
+}
+
+
+
+/* Start with bracket [0,1] with G(0)*G(1)<0.  Search by finding zero crossing
+   sN of line segment [(sL,GL), (sR,GR)].  Set SL=sN if G(sN<0), sR=sN
+   otherwise.*/
+/*---------------------------------------------------------------------------*/
+double
+regulafalsi (double (*G)(double, void*), void *data, double tol, int maxit,
+	     int *iterations)
+/*---------------------------------------------------------------------------*/
+{
+  *iterations = 0;
+
+  double sn = 0.5;//((double *) data)[2]; /* "Undefined" value */
+  double Gn = G(sn, data);
+  if (fabs(Gn) < tol) return sn;
+
+  /* Initial guess is interval [0,1] of course */
+  double s0=0.0;
+  double G0=G(s0, data);
+  if (fabs(G0)<EPS) return s0;
+
+  double s1=1.0;
+  double G1=G(s1, data);
+  if (fabs(G1)<EPS) return s1;
+
+  if (G0*G1 > 0) {
+    print(no_root_str, "regulafalsi", G0, G1);
+    return -1.0;
+  }
+
+  if (G0>G1) {
+    double swap;
+    swap = G0;
+    G0   = G1;
+    G1   = swap;
+  }
+
+  int it = 0;
+  while ( (fabs(Gn) > tol)  &&  (it++ < maxit))  {
+#if 0
+    /* Unmodified Regula-Falsi */
+    /* maintain bracket with G1>G0 */
+    if ( Gn>0 ) {
+      G1 = Gn;
+      s1 = sn;
+    }
+    else {
+      G0 = Gn;
+      s0 = sn;
+    }
+#else
+    /* Modified Regula-Falsi*/
+    if ((Gn>0.0)==(G0>0.0)) {
+      s0 = s1;
+      G0 = G1;
+    }
+    else {
+      // const double gamma_illinois = 0.5;
+      const double gamma_pegasus  = G1/(G1+Gn);
+      G0 *= gamma_pegasus;
+    }
+    s1 = sn;
+    G1 = Gn;
+#endif
+
+    sn = s0 - (s1-s0)*G0/(G1-G0);
+    Gn = G(sn, data);
+  }
+
+  *iterations = it;
+  return sn;
+}
+
+
+
+/* Start with bracket [0,1] with G(0)*G(1)<0.  Search by finding sN=0.5(sL+sR).
+   Set SL=sN if G(sN<0), sR=sN otherwise.*/
+/*---------------------------------------------------------------------------*/
+double
+bisection (double (*G)(double, void*), void *data, double tol, int maxit,
+	   int *iterations)
+/*---------------------------------------------------------------------------*/
+{
+  *iterations = 0;
+
+  double sn = ((double *) data)[2];
+  double Gn = G(sn, data);
+  if (fabs(Gn)<tol) return sn;
+
+  /* Initial guess is interval [0,1] of course */
+  double s0=0.0;
+  double G0=G(s0, data);
+  if (fabs(G0)<EPS) return s0;
+
+  double s1=1.0;
+  double G1=G(s1, data);
+  if (fabs(G1)<EPS) return s1;
+
+  if (G0*G1 > 0.0) {
+    print(no_root_str, "bisection", G0, G1);
+    return -1.0;
+  }
+
+  if (G0>G1) {
+    double swap;
+    swap = G0;
+    G0   = G1;
+    G1   = swap;
+  }
+
+  int it=0;
+  while ( (fabs(Gn)>tol) && (it++ < maxit) ) {
+    if ( Gn>0 ) {
+      G1 = Gn;
+      s1 = sn;
+    }
+    else {
+      G0 = Gn;
+      s0 = sn;
+    }
+
+    sn = 0.5*(s0+s1);
+    Gn = G(sn, data);
+  }
+  *iterations = it;
+  if (it >= maxit) print("Warning: convergence criterion not met\n");
+  return sn;
+}
+
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/nlsolvers.h

@@ -0,0 +1,14 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#ifndef NLSOLVERS_H
+#define NLSOLVERS_H
+
+double bisection   (double (*)(double, void*), void*, double, int, int*);
+double ridder      (double (*)(double, void*), void*, double, int, int*);
+double regulafalsi (double (*)(double, void*), void*, double, int, int*);
+
+#endif /* NLSOLVERS_H */
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/reordersequence.c

@@ -0,0 +1,103 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <opm/core/grid.h>
+#include <opm/core/transport/reorder/reordersequence.h>
+#include <opm/core/transport/reorder/tarjan.h>
+
+
+
+/* Construct adjacency matrix of upwind graph wrt flux.  Column
+   indices are not sorted. */
+static void 
+make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell, 
+                  const double *flux, int *ia, int *ja, int *work)
+{
+    /* Using topology (conn, cptr), and direction, construct adjacency
+       matrix of graph. */
+
+    int i, j, p, f, positive_sign;
+    double theflux;
+
+    /* For each face, store upwind cell in work array */
+    for (i=0; i<nc; ++i) {
+        for (j=faceptr[i]; j<faceptr[i+1]; ++j) {
+            f  = cellfaces[j];
+            positive_sign = (i == face2cell[2*f]);
+            theflux = positive_sign ? flux[f] : -flux[f];
+
+            if ( theflux > 0  ) {
+                /* flux from cell i over face f */
+                work[f] = i; 
+            }
+        }
+    }
+
+    /* Fill ia and ja */
+    p = 0;
+    ia[0] = p;
+    for (i=0; i<nc; ++i) {
+        for (j=faceptr[i]; j<faceptr[i+1]; ++j) {
+
+            f  = cellfaces[j];
+            positive_sign = (i == face2cell[2*f]);
+            theflux = positive_sign ? flux[f] : -flux[f];
+
+            if ( theflux < 0) {
+                ja[p++] = work[f];
+            }
+        }
+
+        ia[i+1] = p;
+    }
+}
+
+static void 
+compute_reorder_sequence(int nc, int nf, int *cellfaces, int *faceptr, int *face2cell, 
+                         const double *flux, int *sequence, int *components, int *ncomponents)
+{
+    int *ia, *ja;
+    int *work;
+    int sz = nf;
+    if (nf < 3*nc) {
+        sz = 3*nc;
+    }
+            
+    work = malloc( sz    * sizeof *work);
+    ia   = malloc((nc+1) * sizeof *ia);
+    ja   = malloc( nf    * sizeof *ja); /* A bit too much... */
+
+    
+    make_upwind_graph(nc, cellfaces, faceptr, face2cell, 
+                      flux, ia, ja, work);
+ 
+    tarjan (nc, ia, ja, sequence, components, ncomponents, work);
+    
+    free(ja);
+    free(ia);
+    free(work);
+}
+
+void compute_sequence(struct UnstructuredGrid *grid, const double *flux, 
+                      int *sequence, 
+                      int *components, int *ncomponents)
+{
+
+    compute_reorder_sequence(grid->number_of_cells, 
+                             grid->number_of_faces, 
+                             grid->cell_faces, 
+                             grid->cell_facepos, 
+                             grid->face_cells,
+                             flux, 
+                             sequence, 
+                             components, 
+                             ncomponents);
+}
+
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */
+
+

opm/core/transport/reorder/reordersequence.h

@@ -0,0 +1,13 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+#ifndef REORDERSEQUENCE_H_INCLUDED
+#define REORDERSEQUENCE_H_INCLUDED
+
+struct UnstructuredGrid;
+void compute_sequence(struct UnstructuredGrid *grid, const double *flux, 
+                      int *sequence, int *components, int *ncomponents);
+
+#endif
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/tarjan.c

@@ -0,0 +1,180 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#include <string.h>
+#include <assert.h>
+#include <opm/core/transport/reorder/tarjan.h>
+
+
+
+
+
+static int min(int a, int b){ return a<b? a : b;}
+
+/* Improved (or obfuscated) version uses less memory
+ *
+ * Use end of P and Q as stack:
+ *   push operation is *s--=elm,
+ *   pop  operation is elm=*++s,
+ *   peek operation is *(s+1)
+ */
+void
+tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
+        int *work)
+/*--------------------------------------------------------------------*/
+{
+    enum {DONE=-2, REMAINING=-1};
+    int c,v,seed,child;
+    int i;
+
+    int *stack    = Q + size,   *bottom  = stack;
+    int *cstack   = P + size-1, *cbottom = cstack;
+
+    int t         = 0;
+    int pos       = 0;
+
+    int   *time   = work;
+    int   *link   = (int *)  time + size;
+    int   *status = (int*)   link + size; /* dual usage... */
+
+  
+
+    memset(work, 0, 3*size    *  sizeof *work);
+    memset(P,    0,   size    *  sizeof *P   );
+    memset(Q,    0,  (size+1) *  sizeof *Q   );
+
+    /* Init status all nodes */
+    for (i=0; i<size; ++i)
+    {
+        status[i] = REMAINING;
+    }
+
+    *ncomp = 0;
+    *Q++ = pos;
+              
+    seed = 0;
+    while (seed < size)
+    { 
+        if (status[seed] == DONE)
+        {
+            ++seed;
+            continue;
+        }
+     
+        *stack-- = seed;     /* push seed */
+
+        t = 0;
+   
+        while ( stack != bottom )
+        {
+            c = *(stack+1);        /* peek c */
+ 
+            assert(status[c] != DONE);
+            assert(status[c] >= -2);
+        
+            if (status[c] == REMAINING)
+            {
+                status[c] = ia[c+1]-ia[c]; /* number of descendants of c */
+                time[c]   = link[c] = t++;
+                *cstack-- = c;             /* push c on strongcomp stack */
+            }
+
+
+
+            /* if all descendants are processed */
+            if (status[c] == 0) 
+            {
+
+                /* if c is root of strong component */
+                if (link[c] == time[c]) 
+                {
+                    do
+                    {
+                        assert (cstack != cbottom);
+
+                        v         = *++cstack; /* pop strong component stack */
+                        status[v] = DONE;
+                        P[pos++]  = v;          /* store vertex in P */
+                    } 
+                    while ( v != c );
+              
+                    *Q++ = pos;       /* store end point of component */ 
+                    ++*ncomp;
+                }
+
+                ++stack;  /* pop c */
+
+                if (stack != bottom) 
+                {
+                    link[*(stack+1)] = min(link[*(stack+1)], link[c]);
+                }
+            }
+
+
+
+            /* if there are more descendants to consider */
+            else 
+            {
+                assert(status[c] > 0);
+
+                child = ja[ia[c] + status[c]-1];
+                --status[c];           /* Decrement descendant count of c*/
+
+                if (status[child] == REMAINING) 
+                {
+                    *stack-- = child; /* push child */
+            
+                } 
+                else if (status[child] >= 0) 
+                {
+                    link[c] = min(link[c], time[child]);
+
+                } 
+                else 
+                {
+                    assert(status[child] = DONE);
+                }
+            }
+        }
+        assert (cstack == cbottom);
+    }
+}
+
+
+#ifdef TEST_TARJAN
+#include <stdio.h>
+#include <stdlib.h>
+/* gcc -DTEST_TARJAN  tarjan.c */
+int main()
+{
+    int size = 5;
+    int ia[] = {0, 1, 2, 4, 5, 5};
+    int ja[] = {1, 2, 0, 3, 4};
+
+    int *work = malloc(3*size*sizeof *work);
+
+    int number_of_components;
+    int *permutation   = malloc(size * sizeof *permutation);
+    int *ptr           = malloc(size * sizeof *ptr);
+    int i,j;
+    tarjan (size, ia, ja, permutation, ptr, &number_of_components, work);
+
+    fprintf(stderr, "Number of components: %d\n", number_of_components);
+    for (i=0; i<number_of_components; ++i)
+    {
+        for (j=ptr[i]; j<ptr[i+1]; ++j)
+        {
+            fprintf(stderr, "%d %d\n", i, permutation[j]);
+        }
+    }
+   
+    free(work);
+    free(permutation);
+    free(ptr);
+   
+    return 0;
+}
+#endif
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/tarjan.h

@@ -0,0 +1,19 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#ifndef TARJAN_H_INCLUDED
+#define TARJAN_H_INCLUDED
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+    void tarjan (int size, int *ia, int *ja, int *rowP, int *P, 
+                 int *ncomp, int *work);
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* TARJAN_H_INCLUDED */
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/twophase.c

@@ -0,0 +1,123 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#include <opm/core/grid.h>
+#include <opm/core/transport/reorder/twophase.h>
+#include <opm/core/transport/reorder/nlsolvers.h>
+
+static struct Parameters get_parameters(struct vdata *vd, 
+                                        const struct cdata *cd, 
+                                        int cell);
+static double fluxfun(double s);
+
+static double 
+fluxfun(double s)
+{
+    return s*s/(s*s + (1-s)*(1-s));
+}
+
+struct Parameters
+{
+    double s0;
+    double influx;  /* sum_j min(v_ij, 0)*f(s_j) */
+    double outflux; /* sum_j max(v_ij, 0)        */
+    double dtpv;    /* dt/pv(i)                  */
+};
+#include <stdio.h>
+static double 
+G(double s, void *data)
+{
+    struct Parameters *p = data;
+    fprintf(stderr, "s:%f dtpv:%f s0:%f, influx:%f outflux:%f G:%f\n", 
+            s, 
+            p->dtpv,
+            p->s0, 
+            p->influx,
+            p->outflux,
+            s - p->s0 +  p->dtpv*(p->outflux*fluxfun(s) + p->influx));
+    /* G(s) = s - s0 + dt/pv*(influx - outflux*f(s) ) */
+    return s - p->s0 +  p->dtpv*(p->outflux*fluxfun(s) + p->influx);
+}
+
+static struct Parameters 
+get_parameters(struct vdata *vd, const struct cdata *cd, int cell)
+{
+    struct UnstructuredGrid *g  = cd->grid;
+    struct Parameters        p;
+
+    p.s0      = vd->saturation[cell];
+    p.influx  = cd->source[cell] >  0 ? -cd->source[cell] : 0.0;
+    p.outflux = cd->source[cell] <= 0 ? -cd->source[cell] : 0.0;
+    p.dtpv    = cd->dt/cd->porevolume[cell];
+    
+    fprintf(stderr, "src:%f pv:%f\n",
+            cd->source[cell],
+            cd->porevolume[cell]);
+
+    int i;
+    vd->saturation[cell] = 0;
+    for (i=g->cell_facepos[cell]; i<g->cell_facepos[cell+1]; ++i) {
+        int    f = g->cell_faces[i];
+        double flux;
+        int    other;
+
+        /* Compute cell flux*/
+        if (cell == g->face_cells[2*f]) {
+            flux  = cd->darcyflux[f];
+            other = g->face_cells[2*f+1];
+        }
+        else {
+            flux  =-cd->darcyflux[f];
+            other = g->face_cells[2*f];
+        }
+       
+        fprintf(stderr, "flux:%f\n", flux);
+        if (other != -1) {
+            if (flux < 0.0) {
+                p.influx  += flux*vd->fractionalflow[other];
+            }
+            else {
+                p.outflux += flux;
+            }
+        }
+    }
+    return p;
+}
+
+
+static enum Method {RIDDER, REGULAFALSI, BISECTION} method = REGULAFALSI;
+
+void solve(void *vdata, const void *cdata, int cell)
+{
+    const struct cdata *cd = cdata;
+    struct vdata       *vd = vdata;
+
+    struct Parameters   p  = get_parameters(vdata, cdata, cell);
+
+    int    iterations;
+    int    maxit = 20;
+    double TOL = 1e-9;
+
+    switch (method) {
+    default:
+    case BISECTION:
+	vd->saturation[cell] = bisection  ( &G, &p, TOL, maxit, &iterations);
+        break;
+
+    case RIDDER:
+	vd->saturation[cell] = ridder     ( &G, &p, TOL, maxit, &iterations); 
+        break;
+
+    case REGULAFALSI:
+	vd->saturation[cell] = regulafalsi( &G, &p, TOL, maxit, &iterations); 
+        break;
+    }
+    
+    vd->fractionalflow[cell]   = fluxfun(vd->saturation[cell]);
+}
+
+
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */
+

opm/core/transport/reorder/twophase.h

@@ -0,0 +1,26 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#ifndef TWOPHASE_H_INCLUDED
+#define TWOPHASE_H_INCLUDED
+
+
+struct vdata {
+    double *saturation;      /* one per cell */
+    double *fractionalflow;  /* one per cell */
+};
+
+struct UnstructuredGrid;
+struct cdata  {
+    struct UnstructuredGrid *grid;
+    const double *darcyflux;   /* one flux per face  in cdata::grid*/
+    const double *source;      /* one source per cell */
+    const double *porevolume;  /* one volume per cell */
+    double  dt;
+};
+
+void solve(void *vdata, const void *cdata, int cell);
+#endif /* TWOPHASE_H_INCLUDED */
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

opm/core/transport/reorder/twophasetransport.c

@@ -0,0 +1,54 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <opm/core/grid.h>
+#include <opm/core/transport/reorder/reordersequence.h>
+#include <opm/core/transport/reorder/tarjan.h>
+#include <opm/core/transport/reorder/twophase.h>
+
+int twophasetransport(
+    const double *porevolume,
+    const double *source,
+    double dt,
+    struct UnstructuredGrid *grid, 
+    const double *darcyflux,
+    double *saturation) 
+{
+    int    sz, i;
+
+    int    ncomponents;
+    int    *sequence;
+    int    *components;
+
+    struct cdata cd = {grid, darcyflux, source, porevolume, dt};
+    struct vdata vd;
+
+    /* Compute sequence of single-cell problems */
+    sequence   = malloc(  grid->number_of_cells    * sizeof *sequence);
+    components = malloc(( grid->number_of_cells+1) * sizeof *components);
+
+    compute_sequence(grid, darcyflux, sequence, components, &ncomponents);
+    assert(ncomponents == grid->number_of_cells);
+
+
+    vd.saturation     = saturation;
+    vd.fractionalflow = malloc(grid->number_of_cells * 
+                               sizeof *vd.fractionalflow);
+    for(i=0; i<grid->number_of_cells; ++i) 
+    {
+        vd.fractionalflow[i] = 0.0;
+    }
+
+    /* Assume all strong components are single-cell domains. */
+    for(i=0; i<grid->number_of_cells; ++i) 
+    { 
+        fprintf(stderr,"i:%d\n", i);
+        solve(&vd, &cd, sequence[i]);
+    }
+
+    free(vd.fractionalflow);
+    free(sequence);
+    free(components);
+}

opm/core/transport/reorder/twophasetransport.h

@@ -0,0 +1,16 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+
+#ifndef TWOPHASETRANSPORT_H_INCLUDED
+#define TWOPHASETRANSPORT_H_INCLUDED
+
+struct UnstructuredGrid;
+int twophasetransport(
+    const double *porevolume,
+    const double *source,
+    double dt,
+    struct UnstructuredGrid *grid, 
+    const double *darcyflux,
+    double *saturation);
+
+#endif /* TWOPHASETRANSPORT_H_INCLUDED */
+