From c0114c52d26ec9452634dbd591d4757a6133dcdb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Thu, 3 Feb 2011 12:51:53 +0100
Subject: [PATCH] Added well handling to cfs_tpfa_expl_mass_transport(),
 simplified interface.

---
 src/cfs_tpfa.c | 53 +++++++++++++++++++++++++++++++++++++++-----------
 src/cfs_tpfa.h | 15 +++++++-------
 2 files changed, 49 insertions(+), 19 deletions(-)

diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c
index a2d05bf9f..94234aaca 100644
--- a/src/cfs_tpfa.c
+++ b/src/cfs_tpfa.c
@@ -1050,25 +1050,37 @@ cfs_tpfa_retrieve_gravtrans(grid_t               *G,
 
 /* ---------------------------------------------------------------------- */
 void
-cfs_tpfa_expl_mass_transport(grid_t       *G,
-                             int           np,
-                             double        dt,
-                             const double *porevol,
-                             const double *masstrans_f,
-                             const double *gravtrans_f,
-                             const double *cpress,
-                             double       *surf_vol)
+cfs_tpfa_expl_mass_transport(grid_t               *G,
+                             well_t               *W,
+                             int                   np,
+                             double                dt,
+                             const double         *porevol,
+                             struct cfs_tpfa_data *h,
+                             double               *surf_vol)
 /* ---------------------------------------------------------------------- */
 {
-    int    c, i, f, c2, p;
-    double dp, dz;
+    int    c, i, f, c2, p, w;
+    double dp, dz, gsgn;
+    const double *masstrans_f, *gravtrans_f, *masstrans_p;
+    const double *cpress, *wpress;
 
+    /* Set up convenience pointers */
+    masstrans_f = h->pimpl->masstrans_f;
+    gravtrans_f = h->pimpl->gravtrans_f;
+    masstrans_p = h->pimpl->masstrans_p;
+    cpress = h->x;
+    wpress = h->x + G->number_of_cells;
+
+    /* Transport through interiour faces */
     for (c = i = 0; c < G->number_of_cells; c++) {
         for (; i < G->cell_facepos[c + 1]; i++) {
             f  = G->cell_faces[i];
 
             if ((c2 = G->face_cells[2*f + 0]) == c) {
+              gsgn = 1.0;
                 c2  = G->face_cells[2*f + 1];
+            } else {
+              gsgn = -1.0;
             }
 
             if (c2 >= 0) {
@@ -1076,13 +1088,32 @@ cfs_tpfa_expl_mass_transport(grid_t       *G,
 
                 for (p = 0; p < np; p++) {
                     dz  = masstrans_f[f*np + p] * dp;
-                    dz += gravtrans_f[f*np + p]     ;
+                    dz += gravtrans_f[f*np + p] * gsgn;
 
+                    /* A positive dz means flow from the cell c into
+                       the cell c2. */
                     surf_vol[c*np + p] -= dz * dt / porevol[c];
                 }
             }
         }
     }
+
+    /* Transport through well perforations */
+    if (W != NULL) {
+        for (w = i = 0; w < W->number_of_wells; w++) {
+            for (; i < W->well_connpos[w + 1]; i++) {
+                c = W->well_cells[i];
+                /* Get pressure difference between cell and well perforation */
+                dp = wpress[i] - cpress[c];
+                for (p = 0; p < np; p++) {
+                    dz = masstrans_p[i*np + p]* dp;
+                    /* A positive dz means flow from the well perforation
+                       i into the cell c */
+                    surf_vol[c*np + p] += dz * dt / porevol[c];
+                }
+            }
+        }
+    }
 }
 
 
diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h
index b9f3997c0..1701c1a37 100644
--- a/src/cfs_tpfa.h
+++ b/src/cfs_tpfa.h
@@ -98,14 +98,13 @@ cfs_tpfa_retrieve_gravtrans(grid_t               *G,
                             double               *gravtrans_f);
 
 void
-cfs_tpfa_expl_mass_transport(grid_t       *G,
-                             int           np,
-                             double        dt,
-                             const double *porevol,
-                             const double *masstrans_f,
-                             const double *gravtrans_f,
-                             const double *cpress,
-                             double       *surf_vol);
+cfs_tpfa_expl_mass_transport(grid_t               *G,
+                             well_t               *W,
+                             int                   np,
+                             double                dt,
+                             const double         *porevol,
+                             struct cfs_tpfa_data *h,
+                             double               *surf_vol);
 
 void
 cfs_tpfa_destroy(struct cfs_tpfa_data *h);