From 634d168b3f759794c617d244d9189c71ff1adfe8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Fri, 25 Nov 2011 14:02:15 +0100
Subject: [PATCH] Assemble well equation form contributions per completion.

We assemble one equation for each well, with the equation being
trivial for wells controlled by bottom-hole targets.  Moreover, at
present we only support BHP and RESV controls.  The infrastructure is
in place to support more general controls (i.e., surface rate
targets), but such wells are currently not part of the 'well.h'
module.
---
 src/cfs_tpfa_residual.c | 91 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 88 insertions(+), 3 deletions(-)

diff --git a/src/cfs_tpfa_residual.c b/src/cfs_tpfa_residual.c
index 47567dc29..c663aa15c 100644
--- a/src/cfs_tpfa_residual.c
+++ b/src/cfs_tpfa_residual.c
@@ -736,6 +736,88 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
 }
 
 
+/* ---------------------------------------------------------------------- */
+static void
+welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
+                 double *res, double *w2c, double *w2w)
+/* ---------------------------------------------------------------------- */
+{
+    int           p;
+    double        fwi;
+    const double *dpflux_w;
+
+    /* Sum reservoir phase flux derivatives set by
+     * compute_darcyflux_and_deriv(). */
+    dpflux_w = h->pimpl->flux_work + (1 * np);
+    for (p = 0, fwi = 0.0; p < np; p++) {
+        fwi += dpflux_w[ p ];
+    }
+
+    *res = fwi * dp;
+    *w2c = 0.0;
+    *w2w = fwi;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static void
+welleq_coeff_resv(int np, struct cfs_tpfa_res_data *h,
+                  double *res, double *w2c, double *w2w)
+/* ---------------------------------------------------------------------- */
+{
+    int           p;
+    const double *pflux, *dpflux_w, *dpflux_c;
+
+     /* Sum reservoir phase flux and its derivatives set by
+      * compute_darcyflux_and_deriv(). */
+    pflux    = h->pimpl->flux_work;
+    dpflux_w = pflux    + (1 * np);
+    dpflux_c = dpflux_w + (1 * np);
+
+    *res = *w2c = *w2w = 0.0;
+    for (p = 0; p < np; p++) {
+        *res += pflux   [ p ];
+        *w2w += dpflux_w[ p ];
+        *w2c += dpflux_c[ p ];
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+static void
+assemble_completion_to_well(int w, int c, int nc, int np,
+                            double pw, double dt,
+                            struct cfs_tpfa_res_wells *W,
+                            struct cfs_tpfa_res_data  *h)
+/* ---------------------------------------------------------------------- */
+{
+    int    wdof;
+    size_t jc, jw;
+    double res, w2c, w2w;
+
+    switch (W->ctrl->ctrl[w]) {
+    case BHP :
+        welleq_coeff_bhp(np, pw - W->ctrl->target[w],
+                         h, &res, &w2c, &w2w);
+        break;
+
+    case RATE:
+        /* Interpret RATE as RESV */
+        welleq_coeff_resv(np, h, &res, &w2c, &w2w);
+        break;
+    }
+
+    /* Assemble completion contributions */
+    wdof = nc + w;
+    jc   = csrmatrix_elm_index(wdof, c   , h->J);
+    jw   = csrmatrix_elm_index(wdof, wdof, h->J);
+
+    h->F    [ wdof ] += dt * res;
+    h->J->sa[ jc   ] += dt * w2c;
+    h->J->sa[ jw   ] += dt * w2w;
+}
+
+
 static void
 assemble_well_contrib(struct cfs_tpfa_res_wells   *W     ,
                       struct compr_quantities_gen *cq    ,
@@ -777,9 +859,12 @@ assemble_well_contrib(struct cfs_tpfa_res_wells   *W     ,
             compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
                                         h->pimpl->flux_work,
                                         h->pimpl->flux_work + np);
-#if 0
-            assemble_completion_to_well();
-#endif
+
+            assemble_completion_to_well(w, c, nc, np, pw, dt, W, h);
+        }
+
+        if (W->ctrl->ctrl[ w ] != BHP) {
+            h->F[ nc + w ] -= dt * W->ctrl->target[ w ];
         }
     }
 }