From 3143333b694e6ccf32bfdc7bcd85d87082c79467 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <thsa@iris.no>
Date: Thu, 8 May 2014 13:32:49 +0200
Subject: [PATCH] Update perfPressure The secondary variable perfpressure is
 updated when the well equation is assembled. The perfpressure is used to
 calculate the welldensity.

---
 opm/autodiff/FullyImplicitBlackoilSolver.hpp      |  1 +
 opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 13 ++++++++++++-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
index 9ec5d7cad..04df36827 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
@@ -170,6 +170,7 @@ namespace Opm {
 
         void
         addWellEq(const SolutionState& state,
+                  WellStateFullyImplicitBlackoil& xw,
                   V& aliveWells);
 
         void updateWellControls(ADB& bhp,
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index ded3d6672..b05706e46 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -712,7 +712,7 @@ namespace {
         // a well control is switched.
         updateWellControls(state.bhp, state.qs, xw);
         V aliveWells;
-        addWellEq(state, aliveWells);
+        addWellEq(state, xw, aliveWells);
         addWellControlEq(state, xw, aliveWells);
     }
 
@@ -722,6 +722,7 @@ namespace {
 
     template <class T>
     void FullyImplicitBlackoilSolver<T>::addWellEq(const SolutionState& state,
+                                                   WellStateFullyImplicitBlackoil& xw,
                                                    V& aliveWells)
     {
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
@@ -904,6 +905,7 @@ namespace {
             residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
         }
 
+
         // Add WELL EQUATIONS
         ADB qs = state.qs;
         for (int phase = 0; phase < np; ++phase) {
@@ -911,6 +913,15 @@ namespace {
 
         }
 
+
+        V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
+        for (int phase = 1; phase < np; ++phase) {
+            cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
+        }
+
+        std::vector<double> cq_d(cq.data(), cq.data() + nperf*np);
+        xw.perfPhaseRates() = cq_d;
+
         residual_.well_flux_eq = qs;
     }