From c3cc4021fa40e4b0d361a42aee4d648f062cd47f Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Mon, 24 Jul 2017 11:31:57 +0200
Subject: [PATCH] adding computeWellRatesWithBhp() to StandardWell

---
 opm/autodiff/StandardWell.hpp      |  6 ++++++
 opm/autodiff/StandardWell_impl.hpp | 32 ++++++++++++++++++++++++++++++
 2 files changed, 38 insertions(+)

diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp
index 356599af0..93c39b117 100644
--- a/opm/autodiff/StandardWell.hpp
+++ b/opm/autodiff/StandardWell.hpp
@@ -244,6 +244,12 @@ namespace Opm
         virtual void wellEqIteration(Simulator& ebosSimulator,
                                      const ModelParameters& param,
                                      WellState& well_state);
+
+        // TODO: maybe we should provide a light version of computeWellFlux, which does not include the
+        // calculation of the derivatives
+        void computeWellRatesWithBhp(const Simulator& ebosSimulator,
+                                     const EvalWell& bhp,
+                                     std::vector<double>& well_flux) const;
     };
 
 }
diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp
index a9973877e..27198db1e 100644
--- a/opm/autodiff/StandardWell_impl.hpp
+++ b/opm/autodiff/StandardWell_impl.hpp
@@ -1936,4 +1936,36 @@ namespace Opm
         updateWellState(xw, param, well_state);
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    computeWellRatesWithBhp(const Simulator& ebosSimulator,
+                            const EvalWell& bhp,
+                            std::vector<double>& well_flux) const
+    {
+        const int np = numberOfPhases();
+        const int numComp = numComponents();
+        well_flux.resize(np, 0.0);
+
+        const bool allow_cf = crossFlowAllowed(ebosSimulator);
+
+        for (int perf = 0; perf < numberOfPerforations(); ++perf) {
+            const int cell_idx = wellCells()[perf];
+            const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
+            // flux for each perforation
+            std::vector<EvalWell> cq_s(numComp, 0.0);
+            std::vector<EvalWell> mob(numComp, 0.0);
+            getMobility(ebosSimulator, perf, mob);
+            computePerfRate(intQuants, mob, wellIndex()[perf], bhp, perfPressureDiffs()[perf], allow_cf, cq_s);
+
+            for(int p = 0; p < np; ++p) {
+                well_flux[p] += cq_s[p].value();
+            }
+        }
+    }
+
 }