From 4b5e1dfadc293bf28613e1aa2806ff36134b44ad Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Tue, 21 Mar 2017 17:04:42 +0100
Subject: [PATCH] more treatment for zero rate wells in
 updateWellStateWithTarget

---
 opm/autodiff/StandardWellsDense_impl.hpp | 86 +++++++++++++++---------
 1 file changed, 56 insertions(+), 30 deletions(-)

diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp
index 20d4b4fb1..46468bcef 100644
--- a/opm/autodiff/StandardWellsDense_impl.hpp
+++ b/opm/autodiff/StandardWellsDense_impl.hpp
@@ -2499,29 +2499,30 @@ namespace Opm {
             // break;
 
         case SURFACE_RATE:
-            // assign target value as initial guess for injectors and
-            // single phase producers (orat, grat, wrat)
+            // checking the number of the phases under control
+            int numPhasesWithTargetsUnderThisControl = 0;
+            for (int phase = 0; phase < np; ++phase) {
+                if (distr[phase] > 0.0) {
+                    numPhasesWithTargetsUnderThisControl += 1;
+                }
+            }
+
+            assert(numPhasesWithTargetsUnderThisControl > 0);
+
             const WellType& well_type = wells().type[well_index];
             if (well_type == INJECTOR) {
+                // assign target value as initial guess for injectors
+                // only handles single phase control at the moment
+                assert(numPhasesWithTargetsUnderThisControl == 1);
+
                 for (int phase = 0; phase < np; ++phase) {
-                    const double& compi = wells().comp_frac[np * well_index + phase];
-                    if (compi > 0.0) {
-                        assert(distr[phase] > 0.);
-                        xw.wellRates()[np*well_index + phase] = target * compi / distr[phase];
+                    if (distr[phase] > 0.) {
+                        xw.wellRates()[np*well_index + phase] = target / distr[phase];
                     } else {
-                        xw.wellRates()[np * well_index + phase] = target * compi;
+                        xw.wellRates()[np * well_index + phase] = 0.;
                     }
                 }
             } else if (well_type == PRODUCER) {
-                // checking the number of the phases under control
-                int numPhasesWithTargetsUnderThisControl = 0;
-                for (int phase = 0; phase < np; ++phase) {
-                    if (distr[phase] > 0.0) {
-                        numPhasesWithTargetsUnderThisControl += 1;
-                    }
-                }
-
-                assert(numPhasesWithTargetsUnderThisControl > 0);
 
                 // update the rates of phases under control based on the target,
                 // and also update rates of phases not under control to keep the rate ratio,
@@ -2540,17 +2541,17 @@ namespace Opm {
                         xw.wellRates()[np * well_index + phase] *= scaling_factor;
                     }
                 } else { // scaling factor is not well defied when orignal_rates_under_phase_control is zero
-                    if (orignal_rates_under_phase_control == 0.0) {
-                        // only handle single-phase control
-                        if (numPhasesWithTargetsUnderThisControl == 1) {
-                            for (int phase = 0; phase < np; ++phase) {
-                                if (distr[phase] > 0.0) {
-                                    xw.wellRates()[np * well_index + phase] = target;
-                                }
-                            }
+                    // separating targets equally between phases under control
+                    const double target_rate_devided = target / numPhasesWithTargetsUnderThisControl;
+                    for (int phase = 0; phase < np; ++phase) {
+                        if (distr[phase] > 0.0) {
+                            xw.wellRates()[np * well_index + phase] = target_rate_devided / distr[phase];
+                        } else {
+                            // this only happens for SURFACE_RATE control
+                            xw.wellRates()[np * well_index + phase] = target_rate_devided;
                         }
                     }
-               }
+                }
             } else {
                 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
             }
@@ -2603,12 +2604,37 @@ namespace Opm {
                 xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
             }
         } else {
-            if (active_[ Water ]) {
-                xw.wellSolutions()[WFrac*nw + well_index] =  wells().comp_frac[np*well_index + Water];
-            }
+            const WellType& well_type = wells().type[well_index];
+            if (well_type == INJECTOR) {
+                // only single phase injection handled
+                if (active_[Water]) {
+                    if (distr[Water] > 0.0) {
+                        xw.wellSolutions()[WFrac * nw + well_index] = 1.0;
+                    } else {
+                        xw.wellSolutions()[WFrac * nw + well_index] = 0.0;
+                    }
+                }
 
-            if (active_[ Gas ]) {
-                xw.wellSolutions()[GFrac*nw + well_index] =  wells().comp_frac[np*well_index + Gas];
+                if (active_[Gas]) {
+                    if (distr[Gas] > 0.0) {
+                        xw.wellSolutions()[GFrac * nw + well_index] = 1.0;
+                    } else {
+                        xw.wellSolutions()[GFrac * nw + well_index] = 0.0;
+                    }
+                }
+
+                // TODO: it is possible to leave injector as a oil well,
+                // when F_w and F_g both equals to zero, not sure under what kind of circumstance
+                // this will happen.
+            } else if (well_type == PRODUCER) { // producers
+                if (active_[Water]) {
+                    xw.wellSolutions()[WFrac * nw + well_index] = 1.0 / np;
+                }
+                if (active_[Gas]) {
+                    xw.wellSolutions()[GFrac * nw + well_index] = 1.0 / np;
+                }
+            } else {
+                OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
             }
         }