From 983a29691c856444baf3b1dde998de1da1bb2cc2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Thu, 23 Aug 2012 14:00:04 +0200
Subject: [PATCH] Function cfs_tpfa_residual_assemble() and friends now return
 singularity flag. The singularity flag is true if there are no pressure
 conditions and no compressibility (so the absolute values of the pressure
 solution will be arbitrary).

---
 opm/core/pressure/tpfa/cfs_tpfa_residual.c | 42 ++++++++++------------
 opm/core/pressure/tpfa/cfs_tpfa_residual.h | 13 +++++--
 2 files changed, 29 insertions(+), 26 deletions(-)

diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c
index fb8ee8795..eb8dc8f34 100644
--- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c
+++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c
@@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid   *G      ,
 
 
 /* ---------------------------------------------------------------------- */
-void
+int
 cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
                       double                       dt       ,
                       struct cfs_tpfa_res_forces  *forces   ,
@@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
                       struct cfs_tpfa_res_data    *h        )
 /* ---------------------------------------------------------------------- */
 {
-    int res_is_neumann, well_is_neumann, c, np2;
+    int res_is_neumann, well_is_neumann, c, np2, singular;
 
     csrmatrix_zero(         h->J);
     vector_zero   (h->J->m, h->F);
@@ -1207,14 +1207,17 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
         assemble_sources(dt, forces->src, h);
     }
 
-    if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) {
-        h->J->sa[0] *= 2;
+    singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp;
+    if (singular) {
+        h->J->sa[0] *= 2.0;
     }
+
+    return singular;
 }
 
 
 /* ---------------------------------------------------------------------- */
-void
+int
 cfs_tpfa_res_comprock_assemble(
                       struct UnstructuredGrid     *G        ,
                       double                       dt       ,
@@ -1240,29 +1243,18 @@ cfs_tpfa_res_comprock_assemble(
      * porevol(pressure)*rock_comp(pressure)/dt
      */
 
-    int     c, w, well_is_neumann, rock_is_incomp;
+    int     c, rock_is_incomp, singular;
     size_t  j;
-    double dpv;
-    const struct Wells* W;
+    double  dpv;
 
     /* Assemble usual system (without rock compressibility). */
-    cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
-                          cpress, wpress, porevol0, h);
-
-    /* Check if we have only Neumann wells. */
-    well_is_neumann = 1;
-    W = forces->wells->W;
-    for (w = 0; well_is_neumann && w < W->number_of_wells; w++) {
-        if ((W->ctrls[w]->current >= 0) && /* OPEN? */
-            (W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) {
-            well_is_neumann = 0;
-        }
-    }
+    singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
+                                     cpress, wpress, porevol0, h);
 
     /* If we made a singularity-removing adjustment in the
        regular assembly, we undo it here. */
-    if (well_is_neumann && h->pimpl->is_incomp) {
-        h->J->sa[0] /= 2;
+    if (singular) {
+        h->J->sa[0] /= 2.0;
     }
 
     /* Add new terms to residual and Jacobian. */
@@ -1280,9 +1272,11 @@ cfs_tpfa_res_comprock_assemble(
     }
 
     /* Re-do the singularity-removing adjustment if necessary */
-    if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) {
-        h->J->sa[0] *= 2;
+    if (rock_is_incomp && singular) {
+        h->J->sa[0] *= 2.0;
     }
+
+    return rock_is_incomp && singular;
 }
 
 
diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h
index 6a86c61bc..5eddeed68 100644
--- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h
+++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h
@@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid   *G      ,
 void
 cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
 
-void
+/* Return value is 1 if the assembled matrix was adjusted to remove a
+   singularity. This happens if all fluids are incompressible and
+   there are no pressure conditions on wells or boundaries.
+   Otherwise return 0.  */
+int
 cfs_tpfa_res_assemble(struct UnstructuredGrid     *G,
                       double                       dt,
                       struct cfs_tpfa_res_forces  *forces,
@@ -72,7 +76,12 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G,
                       const double                *porevol,
                       struct cfs_tpfa_res_data    *h);
 
-void
+/* Return value is 1 if the assembled matrix was adjusted to remove a
+   singularity. This happens if all fluids are incompressible, the
+   rock is incompressible, and there are no pressure conditions on
+   wells or boundaries.
+   Otherwise return 0.  */
+int
 cfs_tpfa_res_comprock_assemble(
                       struct UnstructuredGrid     *G,
                       double                       dt,