Move computation of well potentials from simulator class ot model class

- the computation of well potentials in the model class calculates the
well potentials using computeWellFlux()
- in this way the well potential calculations also handle well where
some perforations are closed by the simulator due to cross-flow.
- the well potentials pr perforation and phase is stored in the well
state.

opm/autodiff/BlackoilModelBase.hpp

@@ -406,6 +406,28 @@ namespace Opm {
                     WellState& well_state);
 
         void
+<<<<<<< HEAD
+=======
+        computeWellFlux(const SolutionState& state,
+                        const std::vector<ADB>& mob_perfcells,
+                        const std::vector<ADB>& b_perfcells,
+                        V& aliveWells,
+                        std::vector<ADB>& cq_s) const;
+
+        void
+        updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
+                                         const SolutionState& state,
+                                         WellState& xw) const;
+
+        void
+        computeWellPotentials(const SolutionState& state,
+                              const std::vector<ADB>& mob_perfcells,
+                              const std::vector<ADB>& b_perfcells,
+                              WellState& well_state);
+
+
+        void
+>>>>>>> Move computation of well potentials from simulator class ot model class
         addWellFluxEq(const std::vector<ADB>& cq_s,
                       const SolutionState& state);
 

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -867,9 +867,104 @@ namespace detail {
         asImpl().addWellFluxEq(cq_s, state);
         asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
         asImpl().addWellControlEq(state, well_state, aliveWells);
+
+        asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state);
+
     }
 
+    template <class Grid, class Implementation>
+    void
+    BlackoilModelBase<Grid, Implementation>::
+    computeWellPotentials(const SolutionState& state,
+                          const std::vector<ADB>& mob_perfcells,
+                          const std::vector<ADB>& b_perfcells,
+                          WellState& well_state)
+    {
+        const int nw = wells().number_of_wells;
+        const int np = wells().number_of_phases;
+//        const Opm::PhaseUsage pu = fluid_.phaseUsage();
+        V bhps = V::Zero(nw);
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* ctrl = wells().ctrls[w];
+            const int nwc = well_controls_get_num(ctrl);
+            //Loop over all controls until we find a BHP control
+            //or a THP control that specifies what we need.
+            //Pick the value that gives largest potential flow
+            for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
 
+                if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
+                    bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
+                }
+
+//                if(well_controls_iget_type(ctrl, ctrl_index) == THP) {
+//                    double aqua = 0.0;
+//                    double liquid = 0.0;
+//                    double vapour = 0.0;
+
+//                    if (active_[ Water ]) {
+//                        aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
+//                    }
+//                    if (active_[ Oil ]) {
+//                        liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
+//                    }
+//                    if (active_[ Gas ]) {
+//                        vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
+//                    }
+
+//                    const int vfp        = well_controls_iget_vfp(ctrl, ctrl_index);
+//                    const double& thp    = well_controls_iget_target(ctrl, ctrl_index);
+//                    const double& alq    = well_controls_iget_alq(ctrl, ctrl_index);
+
+//                    //Set *BHP* target by calculating bhp from THP
+//                    const WellType& well_type = wells().type[w];
+
+//                    if (well_type == INJECTOR) {
+//                        double dp = detail::computeHydrostaticCorrection(
+//                                    wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
+//                                    well_perforation_densities_, gravity);
+//                        const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
+//                        // pick the bhp that gives the largest potentials i.e. largest bhp for injectors
+//                        if ( bhp > bhps[w]) {
+//                            bhps[w] = bhp;
+//                        }
+//                    }
+//                    else if (well_type == PRODUCER) {
+//                        double dp = detail::computeHydrostaticCorrection(
+//                                    wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
+//                                    well_perforation_densities_, gravity);
+
+//                        const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
+//                        // pick the bhp that gives the largest potentials i.e. smalest bhp for producers
+//                        if ( bhp < bhps[w]) {
+//                            bhps[w] = bhp;
+//                        }
+//                    }
+//                    else {
+//                        OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
+//                    }
+//                }
+            }
+        }
+
+        // use bhp limit from control
+        SolutionState state0 = state;
+        asImpl().makeConstantState(state0);
+        state0.bhp = ADB::constant(bhps);
+
+        // compute well potentials
+        V aliveWells;
+        std::vector<ADB> well_potentials;
+        asImpl().computeWellFlux(state0, mob_perfcells,  b_perfcells, aliveWells, well_potentials);
+
+        // store well potentials in the well state
+        // transform to a single vector instead of separate vectors pr phase
+        const int nperf = wells().well_connpos[nw];
+        V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
+        for (int phase = 1; phase < np; ++phase) {
+            cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
+        }
+        well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
+    }
 
 
 

opm/autodiff/SimulatorBase_impl.hpp

@@ -396,30 +396,12 @@ namespace Opm
     {
         const int nw = wells->number_of_wells;
         const int np = wells->number_of_phases;
-
         well_potentials.clear();
-        well_potentials.resize(nw*np,0.0);       
+        well_potentials.resize(nw*np,0.0);
         for (int w = 0; w < nw; ++w) {
             for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
-                const double well_cell_pressure = x.pressure()[wells->well_cells[perf]];
-                const double drawdown_used = well_cell_pressure - xw.perfPress()[perf];
-                const WellControls* ctrl = wells->ctrls[w];
-                const int nwc = well_controls_get_num(ctrl);
-                //Loop over all controls until we find a BHP control
-                //that specifies what we need...
-                double bhp = 0.0;
-                for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
-                    if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
-                        bhp = well_controls_iget_target(ctrl, ctrl_index);
-                    }
-                    // TODO: do something for thp;
-                }
-                // Calculate the pressure difference in the well perforation
-                const double dp = xw.perfPress()[perf] - xw.bhp()[w];
-                const double drawdown_maximum = well_cell_pressure - (bhp + dp);
-
                 for (int phase = 0; phase < np; ++phase) {
-                    well_potentials[w*np + phase] += (drawdown_maximum / drawdown_used * xw.perfPhaseRates()[perf*np + phase]);
+                    well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
                 }
             }
         }

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -102,6 +102,9 @@ namespace Opm
                 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
             }
 
+            well_potentials_.clear();
+            well_potentials_.resize(nperf * np, 0.0);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if( ! prevState.wellMap().empty() )
@@ -184,9 +187,14 @@ namespace Opm
         std::vector<int>& currentControls() { return current_controls_; }
         const std::vector<int>& currentControls() const { return current_controls_; }
 
+        /// One rate per phase and well connection.
+        std::vector<double>& wellPotentials() { return well_potentials_; }
+        const std::vector<double>& wellPotentials() const { return well_potentials_; }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
+        std::vector<double> well_potentials_;
     };
 
 } // namespace Opm