Merge remote-tracking branch 'origin/master' into frankenstein

opm/autodiff/BlackoilModelBase.hpp

@@ -32,6 +32,7 @@
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/autodiff/BlackoilModelEnums.hpp>
 #include <opm/autodiff/VFPProperties.hpp>
+#include <opm/autodiff/RateConverter.hpp>
 #include <opm/autodiff/IterationReport.hpp>
 #include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
@@ -113,6 +114,10 @@ namespace Opm {
         typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
         typedef typename ModelTraits<Implementation>::SolutionState SolutionState;
 
+        // For the conversion between the surface volume rate and resrevoir voidage rate
+        using RateConverterType = RateConverter::
+                                  SurfaceToReservoirVoidage<BlackoilPropsAdInterface, std::vector<int> >;
+
         // ---------  Public methods  ---------
 
         /// Construct the model. It will retain references to the
@@ -267,6 +272,18 @@ namespace Opm {
         computeFluidInPlace(const ReservoirState& x,
                             const std::vector<int>& fipnum);
 
+        /// Function to compute the resevoir voidage for the production wells.
+        /// TODO: Probably should go to well model, while we then have duplications there for two Well Models.
+        /// With time, it looks like probably we will introduce a base class for Well Models.
+        void computeWellVoidageRates(const ReservoirState& reservoir_state,
+                                     const WellState& well_state,
+                                     std::vector<double>& well_voidage_rates,
+                                     std::vector<double>& voidage_conversion_coeffs);
+
+
+        void applyVREPGroupControl(const ReservoirState& reservoir_state,
+                                   WellState& well_state);
+
     protected:
 
         // ---------  Types and enums  ---------
@@ -321,6 +338,9 @@ namespace Opm {
         double current_relaxation_;
         V dx_old_;
 
+        // rate converter between the surface volume rates and reservoir voidage rates
+        RateConverterType rate_converter_;
+
         // ---------  Protected methods  ---------
 
         /// Access the most-derived class used for
@@ -379,6 +399,7 @@ namespace Opm {
         IterationReport
         solveWellEq(const std::vector<ADB>& mob_perfcells,
                     const std::vector<ADB>& b_perfcells,
+                    const ReservoirState& reservoir_state,
                     SolutionState& state,
                     WellState& well_state);
 

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -137,6 +137,11 @@ typedef Eigen::Array<double,
         , terminal_output_ (terminal_output)
         , material_name_(0)
         , current_relaxation_(1.0)
+        // only one region 0 used, which means average reservoir hydrocarbon conditions in
+        // the field will be calculated.
+        // TODO: more delicate implementation will be required if we want to handle different
+        // FIP regions specified from the well specifications.
+        , rate_converter_(fluid_, std::vector<int>(AutoDiffGrid::numCells(grid_),0))
     {
         if (active_[Water]) {
             material_name_.push_back("Water");
@@ -719,6 +724,13 @@ typedef Eigen::Array<double,
         // get reasonable initial conditions for the wells
         asImpl().wellModel().updateWellControls(well_state);
 
+        if (asImpl().wellModel().wellCollection()->groupControlActive()) {
+            // enforce VREP control when necessary.
+            applyVREPGroupControl(reservoir_state, well_state);
+
+            asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
+        }
+
         // Create the primary variables.
         SolutionState state = asImpl().variableState(reservoir_state, well_state);
 
@@ -755,7 +767,7 @@ typedef Eigen::Array<double,
         asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
-            iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+            iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
         }
         V aliveWells;
         std::vector<ADB> cq_s;
@@ -770,6 +782,7 @@ typedef Eigen::Array<double,
             asImpl().makeConstantState(state0);
             asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
         }
+
         return iter_report;
     }
 
@@ -901,8 +914,10 @@ typedef Eigen::Array<double,
         // Add well contributions to mass balance equations
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
         const int np = asImpl().numPhases();
+        const V& efficiency_factors = wellModel().wellPerfEfficiencyFactors();
         for (int phase = 0; phase < np; ++phase) {
-            residual_.material_balance_eq[phase] -= superset(cq_s[phase], wellModel().wellOps().well_cells, nc);
+            residual_.material_balance_eq[phase] -= superset(efficiency_factors * cq_s[phase],
+                                                             wellModel().wellOps().well_cells, nc);
         }
     }
 
@@ -951,6 +966,7 @@ typedef Eigen::Array<double,
     BlackoilModelBase<Grid, WellModel, Implementation>::
     solveWellEq(const std::vector<ADB>& mob_perfcells,
                 const std::vector<ADB>& b_perfcells,
+                const ReservoirState& reservoir_state,
                 SolutionState& state,
                 WellState& well_state)
     {
@@ -1017,6 +1033,12 @@ typedef Eigen::Array<double,
             // wells active or not as parallel logging will take place that needs to
             // communicate with all processes.
             asImpl().wellModel().updateWellControls(well_state);
+
+            if (asImpl().wellModel().wellCollection()->groupControlActive()) {
+                // Enforce the VREP control
+                applyVREPGroupControl(reservoir_state, well_state);
+                asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
+            }
         } while (it < 15);
 
         if (converged) {
@@ -1700,6 +1722,7 @@ typedef Eigen::Array<double,
         const double residualWell     = detail::infinityNormWell(residual_.well_eq,
                                                                  linsolver_.parallelInformation());
         converged_Well = converged_Well && (residualWell < tol_well_control);
+
         const bool converged = converged_MB && converged_CNV && converged_Well;
 
         // Residual in Pascal can have high values and still be ok.
@@ -1820,6 +1843,7 @@ typedef Eigen::Array<double,
         const double residualWell     = detail::infinityNormWell(residual_.well_eq,
                                                                  linsolver_.parallelInformation());
         converged_Well = converged_Well && (residualWell < tol_well_control);
+
         const bool converged = converged_Well;
 
         // if one of the residuals is NaN, throw exception, so that the solver can be restarted
@@ -2360,6 +2384,117 @@ typedef Eigen::Array<double,
         return values;
     }
 
+
+
+
+
+    template <class Grid, class WellModel, class Implementation>
+    void
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeWellVoidageRates(const ReservoirState& reservoir_state,
+                            const WellState& well_state,
+                            std::vector<double>& well_voidage_rates,
+                            std::vector<double>& voidage_conversion_coeffs)
+    {
+        // TODO: for now, we store the voidage rates for all the production wells.
+        // For injection wells, the rates are stored as zero.
+        // We only store the conversion coefficients for all the injection wells.
+        // Later, more delicate model will be implemented here.
+        // And for the moment, group control can only work for serial running.
+        const int nw = well_state.numWells();
+        const int np = numPhases();
+
+        const Wells* wells = asImpl().wellModel().wellsPointer();
+
+        // we calculate the voidage rate for each well, that means the sum of all the phases.
+        well_voidage_rates.resize(nw, 0);
+        // store the conversion coefficients, while only for the use of injection wells.
+        voidage_conversion_coeffs.resize(nw * np, 1.0);
+
+        int global_number_wells = nw;
+
+#if HAVE_MPI
+        if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+        {
+            const auto& info =
+                boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+            global_number_wells = info.communicator().sum(global_number_wells);
+            if ( global_number_wells )
+            {
+                rate_converter_.defineState(reservoir_state, boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()));
+            }
+        }
+        else
+#endif
+        {
+            if ( global_number_wells )
+            {
+                rate_converter_.defineState(reservoir_state);
+            }
+        }
+
+        std::vector<double> well_rates(np, 0.0);
+        std::vector<double> convert_coeff(np, 1.0);
+
+
+        if ( !well_voidage_rates.empty() ) {
+            for (int w = 0; w < nw; ++w) {
+                const bool is_producer = wells->type[w] == PRODUCER;
+
+                // not sure necessary to change all the value to be positive
+                if (is_producer) {
+                    std::transform(well_state.wellRates().begin() + np * w,
+                                   well_state.wellRates().begin() + np * (w + 1),
+                                   well_rates.begin(), std::negate<double>());
+
+                    // the average hydrocarbon conditions of the whole field will be used
+                    const int fipreg = 0; // Not considering FIP for the moment.
+
+                    rate_converter_.calcCoeff(well_rates, fipreg, convert_coeff);
+                    well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
+                                                               convert_coeff.begin(), 0.0);
+                } else {
+                    // TODO: Not sure whether will encounter situation with all zero rates
+                    // and whether it will cause problem here.
+                    std::copy(well_state.wellRates().begin() + np * w,
+                              well_state.wellRates().begin() + np * (w + 1),
+                              well_rates.begin());
+                    // the average hydrocarbon conditions of the whole field will be used
+                    const int fipreg = 0; // Not considering FIP for the moment.
+                    rate_converter_.calcCoeff(well_rates, fipreg, convert_coeff);
+                    std::copy(convert_coeff.begin(), convert_coeff.end(),
+                              voidage_conversion_coeffs.begin() + np * w);
+                }
+            }
+        }
+    }
+
+
+
+
+
+    template <class Grid, class WellModel, class Implementation>
+    void
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    applyVREPGroupControl(const ReservoirState& reservoir_state,
+                          WellState& well_state)
+    {
+        if (asImpl().wellModel().wellCollection()->havingVREPGroups()) {
+            std::vector<double> well_voidage_rates;
+            std::vector<double> voidage_conversion_coeffs;
+            computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs);
+            asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
+
+            // for the wells under group control, update the currentControls for the well_state
+            for (const WellNode* well_node : asImpl().wellModel().wellCollection()->getLeafNodes()) {
+                if (well_node->isInjector() && !well_node->individualControl()) {
+                    const int well_index = well_node->selfIndex();
+                    well_state.currentControls()[well_index] = well_node->groupControlIndex();
+                }
+            }
+        }
+    }
+
 } // namespace Opm
 
 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -169,6 +169,7 @@ namespace Opm {
         IterationReport
         solveWellEq(const std::vector<ADB>& mob_perfcells,
                     const std::vector<ADB>& b_perfcells,
+                    const ReservoirState& reservoir_state,
                     SolutionState& state,
                     WellState& well_state);
 

opm/autodiff/BlackoilMultiSegmentModel_impl.hpp

@@ -150,6 +150,11 @@ namespace Opm {
         // get reasonable initial conditions for the wells
         wellModel().updateWellControls(well_state);
 
+        // enforce VREP control when necessary.
+        Base::applyVREPGroupControl(reservoir_state, well_state);
+
+        asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
+
         // Create the primary variables.
         SolutionState state = asImpl().variableState(reservoir_state, well_state);
 
@@ -198,7 +203,7 @@ namespace Opm {
         wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
-            iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+            iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
         }
 
         // the perforation flux here are different
@@ -221,10 +226,11 @@ namespace Opm {
     IterationReport
     BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
                                                  const std::vector<ADB>& b_perfcells,
+                                                 const ReservoirState& reservoir_state,
                                                  SolutionState& state,
                                                  WellState& well_state)
     {
-        IterationReport iter_report = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+        IterationReport iter_report = Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
 
         if (iter_report.converged) {
             // We must now update the state.segp and state.segqs members,

opm/autodiff/BlackoilTransportModel.hpp

@@ -135,7 +135,7 @@ namespace Opm {
             asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
             if (param_.solve_welleq_initially_ && initial_assembly) {
                 // solve the well equations as a pre-processing step
-                iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+                iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
             }
             V aliveWells;
             std::vector<ADB> cq_s;

opm/autodiff/MultisegmentWells.cpp

@@ -141,10 +141,13 @@ namespace Opm {
 
     MultisegmentWells::
     MultisegmentWells(const Wells* wells_arg,
+                      WellCollection* well_collection,
                       const std::vector< const Well* >& wells_ecl,
                       const int time_step)
       : wells_multisegment_( createMSWellVector(wells_arg, wells_ecl, time_step) )
       , wops_ms_(wells_multisegment_)
+      , well_collection_(well_collection)
+      , well_perforation_efficiency_factors_(Vector::Ones(numWells()))
       , num_phases_(wells_arg ? wells_arg->number_of_phases : 0)
       , wells_(wells_arg)
       , fluid_(nullptr)
@@ -270,6 +273,8 @@ namespace Opm {
         }
 
         assert(start_perforation == nperf_total_);
+
+        calculateEfficiencyFactors();
     }
 
 
@@ -381,6 +386,56 @@ namespace Opm {
         return indices;
     }
 
+
+
+
+
+    WellCollection*
+    MultisegmentWells::
+    wellCollection() const {
+        return well_collection_;
+    }
+
+
+
+
+
+    void
+    MultisegmentWells::
+    calculateEfficiencyFactors()
+    {
+        if ( !localWellsActive() ) {
+            return;
+        }
+        // get efficiency factor for each well first
+        const int nw = wells_->number_of_wells;
+
+        Vector well_efficiency_factors = Vector::Ones(nw);
+
+        for (int w = 0; w < nw; ++w) {
+            const std::string well_name = wells_->name[w];
+            // get the well node in the well collection
+            WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
+            well_efficiency_factors(w) = well_node.getAccumulativeEfficiencyFactor();
+        }
+
+        // map them to the perforation.
+        well_perforation_efficiency_factors_ = wellOps().w2p * well_efficiency_factors.matrix();
+    }
+
+
+
+
+
+    const
+    MultisegmentWells::Vector&
+    MultisegmentWells::
+    wellPerfEfficiencyFactors() const
+    {
+        return well_perforation_efficiency_factors_;
+    }
+
+
 } // end of namespace Opm
 
 

opm/autodiff/MultisegmentWells.hpp

@@ -32,6 +32,7 @@
 #include <cassert>
 
 #include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/wells/WellCollection.hpp>
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -89,6 +90,7 @@ namespace Opm {
             // TODO: using a vector of WellMultiSegmentConstPtr for now
             // TODO: it should use const Wells or something else later.
             MultisegmentWells(const Wells* wells_arg,
+                              WellCollection* well_collection,
                               const std::vector< const Well* >& wells_ecl,
                               const int time_step);
 
@@ -229,11 +231,26 @@ namespace Opm {
                                                 const std::vector<ADB>& kr_adb,
                                                 const std::vector<ADB>& fluid_density);
 
+            WellCollection* wellCollection() const;
+
+            void calculateEfficiencyFactors();
+
+            const Vector& wellPerfEfficiencyFactors() const;
+
+
     protected:
         // TODO: probably a wells_active_ will be required here.
         bool wells_active_;
         std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
         MultisegmentWellOps wops_ms_;
+        // It will probably need to be updated during running time.
+        WellCollection* well_collection_;
+
+        // The efficiency factor for each connection
+        // It is specified based on wells and groups
+        // By default, they should all be one.
+        Vector well_perforation_efficiency_factors_;
+
         const int num_phases_;
         int nseg_total_;
         int nperf_total_;

opm/autodiff/MultisegmentWells_impl.hpp

@@ -904,6 +904,20 @@ namespace Opm
                 break;
             }
 
+            if (wellCollection()->groupControlActive()) {
+                // get the well node in the well collection
+                WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
+
+                // update whehter the well is under group control or individual control
+                if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
+                    // under group control
+                    well_node.setIndividualControl(false);
+                } else {
+                    // individual control
+                    well_node.setIndividualControl(true);
+                }
+            }
+
         }
 
     }

opm/autodiff/SimulatorBase_impl.hpp

@@ -201,7 +201,7 @@ namespace Opm
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
 
-            const WellModel well_model(wells);
+            const WellModel well_model(wells, &(wells_manager.wellCollection()));
 
             std::unique_ptr<Solver> solver = asImpl().createSolver(well_model);
 

opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp

@@ -122,7 +122,7 @@ namespace Opm
             const auto wells_ecl = eclipse_state_->getSchedule().getWells(timer.currentStepNum());
             const int current_time_step = timer.currentStepNum();
 
-            const WellModel well_model(wells, wells_ecl, current_time_step);
+            const WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, current_time_step);
 
             well_state.init(well_model, state, prev_well_state, wells);
 

opm/autodiff/StandardWells.hpp

@@ -38,6 +38,7 @@
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/DynamicListEconLimited.hpp>
+#include <opm/core/wells/WellCollection.hpp>
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
@@ -70,7 +71,7 @@ namespace Opm {
                                             Eigen::Dynamic,
                                             Eigen::RowMajor>;
             // ---------  Public methods  ---------
-            explicit StandardWells(const Wells* wells_arg);
+            StandardWells(const Wells* wells_arg, WellCollection* well_collection);
 
             void init(const BlackoilPropsAdInterface* fluid_arg,
                       const std::vector<bool>* active_arg,
@@ -190,10 +191,24 @@ namespace Opm {
                                   const WellState& well_state,
                                   DynamicListEconLimited& list_econ_limited) const;
 
+
+            WellCollection* wellCollection() const;
+
+            void calculateEfficiencyFactors();
+
+            const Vector& wellPerfEfficiencyFactors() const;
+
         protected:
             bool wells_active_;
             const Wells*   wells_;
             const WellOps  wops_;
+            // It will probably need to be updated during running time.
+            WellCollection* well_collection_;
+
+            // The efficiency factor for each connection
+            // It is specified based on wells and groups
+            // By default, they should all be one.
+            Vector well_perforation_efficiency_factors_;
 
             const BlackoilPropsAdInterface* fluid_;
             const std::vector<bool>*  active_;

opm/autodiff/StandardWellsSolvent.hpp

@@ -37,7 +37,7 @@ namespace Opm {
             using Base::computeWellConnectionDensitesPressures;
 
             // ---------  Public methods  ---------
-            StandardWellsSolvent(const Wells* wells_arg);
+            StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection);
 
             // added the Solvent related
             void initSolvent(const SolventPropsAdFromDeck* solvent_props,

opm/autodiff/StandardWellsSolvent_impl.hpp

@@ -32,8 +32,8 @@ namespace Opm
 
 
 
-    StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg)
+    StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection)
-        : Base(wells_arg)
+        : Base(wells_arg, well_collection)
         , solvent_props_(nullptr)
         , solvent_pos_(-1)
         , has_solvent_(false)

opm/autodiff/StandardWells_impl.hpp

@@ -71,10 +71,12 @@ namespace Opm
 
 
 
-StandardWells::StandardWells(const Wells* wells_arg)
+    StandardWells::StandardWells(const Wells* wells_arg, WellCollection* well_collection)
       : wells_active_(wells_arg!=nullptr)
       , wells_(wells_arg)
       , wops_(wells_arg)
+      , well_collection_(well_collection)
+      , well_perforation_efficiency_factors_(Vector::Ones(wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0))
       , fluid_(nullptr)
       , active_(nullptr)
       , phase_condition_(nullptr)
@@ -103,6 +105,8 @@ StandardWells::StandardWells(const Wells* wells_arg)
         vfp_properties_ = vfp_properties_arg;
         gravity_ = gravity_arg;
         perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
+
+        calculateEfficiencyFactors();
     }
 
 
@@ -741,6 +745,7 @@ StandardWells::StandardWells(const Wells* wells_arg)
                     break;
                 }
             }
+
             if (ctrl_index != nwc) {
                 // Constraint number ctrl_index was broken, switch to it.
                 // We disregard terminal_ouput here as with it only messages
@@ -748,6 +753,7 @@ StandardWells::StandardWells(const Wells* wells_arg)
                 logger.wellSwitched(wells().name[w],
                                     well_controls_iget_type(wc, current),
                                     well_controls_iget_type(wc, ctrl_index));
+
                 xw.currentControls()[w] = ctrl_index;
                 current = xw.currentControls()[w];
             }
@@ -847,6 +853,22 @@ StandardWells::StandardWells(const Wells* wells_arg)
 
                 break;
             }
+
+
+            if (wellCollection()->groupControlActive()) {
+
+                // get well node in the well collection
+                WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
+
+                // update whehter the well is under group control or individual control
+                if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
+                    // under group control
+                    well_node.setIndividualControl(false);
+                } else {
+                    // individual control
+                    well_node.setIndividualControl(true);
+                }
+            }
         }
 
     }
@@ -1547,4 +1569,48 @@ StandardWells::StandardWells(const Wells* wells_arg)
     }
 
 
+
+
+
+    WellCollection* StandardWells::wellCollection() const
+    {
+        return well_collection_;
+    }
+
+
+
+
+
+    void StandardWells::calculateEfficiencyFactors()
+    {
+        if ( !localWellsActive() ) {
+            return;
+        }
+        // get efficiency factor for each well first
+        const int nw = wells_->number_of_wells;
+
+        Vector well_efficiency_factors = Vector::Ones(nw);
+
+        for (int w = 0; w < nw; ++w) {
+            const std::string well_name = wells_->name[w];
+            const WellNode& well_node = well_collection_->findWellNode(well_name);
+
+            well_efficiency_factors(w) = well_node.getAccumulativeEfficiencyFactor();
+        }
+
+        // map them to the perforation.
+        well_perforation_efficiency_factors_ = wellOps().w2p * well_efficiency_factors.matrix();
+    }
+
+
+
+
+
+    const StandardWells::Vector&
+    StandardWells::wellPerfEfficiencyFactors() const
+    {
+        return well_perforation_efficiency_factors_;
+    }
+
+
 } // namespace Opm

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -166,8 +166,10 @@ namespace Opm
                             // If the set of controls have changed, this may not be identical
                             // to the last control, but it must be a valid control.
                             currentControls()[ newIndex ] = old_control_index;
+                        } else {
+                            assert(well_controls_get_num(wells->ctrls[w]) > 0);
+                            currentControls()[ newIndex ] = 0;
                         }
-
                     }
                 }
             }

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -542,7 +542,7 @@ namespace Opm {
         wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
-            Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+            Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
         }
 
         V aliveWells;

tests/test_multisegmentwells.cpp

@@ -119,7 +119,7 @@ struct SetupMSW {
         const Wells* wells = wells_manager.c_wells();
         const auto wells_ecl = ecl_state.getSchedule().getWells(current_timestep);
 
-        ms_wells.reset(new Opm::MultisegmentWells(wells, wells_ecl, current_timestep));
+        ms_wells.reset(new Opm::MultisegmentWells(wells, &(wells_manager.wellCollection()), wells_ecl, current_timestep));
     };
 
     std::shared_ptr<const Opm::MultisegmentWells> ms_wells;