From b83f37dcc0530df36bbaee9e33ccfca4e3893b8a Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Thu, 23 Mar 2017 16:36:48 +0100
Subject: [PATCH] considering the distr has zero values for non-injecting
 phases

to avoid NaN resulting from devided by zero.
---
 opm/autodiff/StandardWellsDense_impl.hpp | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp
index 37fc38d24..c8f2efb27 100644
--- a/opm/autodiff/StandardWellsDense_impl.hpp
+++ b/opm/autodiff/StandardWellsDense_impl.hpp
@@ -1282,7 +1282,11 @@ namespace Opm {
             if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
                 const double* distr = well_controls_iget_distr(wc, current);
                 for (int p = 0; p < np; ++p) {
-                    F[p] /= distr[p];
+                    if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value
+                        F[p] /= distr[p];
+                    } else {
+                        F[p] = 0.;
+                    }
                 }
             } else {
                 for (int p = 0; p < np; ++p) {
@@ -2201,7 +2205,13 @@ namespace Opm {
         const WellControls* wc = wells().ctrls[wellIdx];
         if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
             const double* distr = well_controls_get_current_distr(wc);
-            return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
+            if (distr[phaseIdx] > 0.) {
+                return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
+            } else {
+                // TODO: not sure why return EvalWell(0.) causing problem here
+                // Probably due to the wrong Jacobians.
+                return wellVolumeFraction(wellIdx, phaseIdx);
+            }
         }
         std::vector<double> g = {1,1,0.01};
         return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]);