revert changes

opm/simulators/wells/MultisegmentWell.hpp

@@ -132,6 +132,8 @@ namespace Opm
                                                  const WellState& well_state,
                                                  DeferredLogger& deferred_logger) override; // should be const?
 
+        void updateIPRImplicit(const Simulator& ebos_simulator, WellState& well_state, DeferredLogger& deferred_logger);                                                 
+
         virtual void updateProductivityIndex(const Simulator& ebosSimulator,
                                              const WellProdIndexCalculator& wellPICalc,
                                              WellState& well_state,
@@ -246,6 +248,10 @@ namespace Opm
                                  const Simulator& ebos_simulator,
                                  DeferredLogger& deferred_logger) const;
 
+        bool computeWellPotentialsImplicit(const Simulator& ebos_simulator,
+                                           std::vector<double>& well_potentials,
+                                           DeferredLogger& deferred_logger) const;                                      
+
         virtual double getRefDensity() const override;
 
         virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
@@ -262,7 +268,8 @@ namespace Opm
                                                 const Well::ProductionControls& prod_controls,
                                                 WellState& well_state,
                                                 const GroupState& group_state,
-                                                DeferredLogger& deferred_logger) override;
+                                                DeferredLogger& deferred_logger, 
+                                                const bool allow_switch = true) override;
 
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                                     const double dt,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -286,13 +286,23 @@ namespace Opm
         }
 
         debug_cost_counter_ = 0;
-        // does the well have a THP related constraint?
+        bool converged_implicit = false;
-        const auto& summaryState = ebosSimulator.vanguard().summaryState();
+        if (this->param_.local_well_solver_control_switching_) {
-        if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
+            converged_implicit = computeWellPotentialsImplicit(ebosSimulator, well_potentials, deferred_logger);
-            computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
+            if (!converged_implicit) {
-        } else {
+                deferred_logger.debug("Implicit potential calculations failed for well "
-            well_potentials = computeWellPotentialWithTHP(
+                                       + this->name() + ",  reverting to original aproach.");
-                well_state, ebosSimulator, deferred_logger);
+            }
+        }
+        if (!converged_implicit) {
+            // does the well have a THP related constraint?
+            const auto& summaryState = ebosSimulator.vanguard().summaryState();
+            if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
+                computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
+            } else {
+                well_potentials = computeWellPotentialWithTHP(
+                    well_state, ebosSimulator, deferred_logger);
+            }
         }
         deferred_logger.debug("Cost in iterations of finding well potential for well "
                               + this->name() + ": " + std::to_string(debug_cost_counter_));
@@ -488,7 +498,70 @@ namespace Opm
         return potentials;
     }
 
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    computeWellPotentialsImplicit(const Simulator& ebos_simulator,
+                                  std::vector<double>& well_potentials,
+                                  DeferredLogger& deferred_logger) const
+    {
+        // Create a copy of the well.
+        // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials 
+        // is allready a copy, but not from other calls. 
+        MultisegmentWell<TypeTag> well_copy(*this);
+        well_copy.debug_cost_counter_ = 0;
 
+        // store a copy of the well state, we don't want to update the real well state
+        WellState well_state_copy = ebos_simulator.problem().wellModel().wellState();
+        const auto& group_state = ebos_simulator.problem().wellModel().groupState();
+        auto& ws = well_state_copy.well(this->index_of_well_);
+        
+        // get current controls        
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        auto inj_controls = well_copy.well_ecl_.isInjector()
+            ? well_copy.well_ecl_.injectionControls(summary_state)
+            : Well::InjectionControls(0);
+        auto prod_controls = well_copy.well_ecl_.isProducer()
+            ? well_copy.well_ecl_.productionControls(summary_state)
+            : Well::ProductionControls(0);
+        
+        // prepare/modify well state and control
+        well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
+        
+        well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
+
+        // initialize rates from previous potentials
+        const int np = this->number_of_phases_;
+        bool trivial = true;
+        for (int phase = 0; phase < np; ++phase){
+            trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
+        }
+        if (!trivial) {
+            const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
+            for (int phase = 0; phase < np; ++phase) {
+                ws.surface_rates[phase] = sign * ws.well_potentials[phase];
+            }
+        }
+        well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
+                                                 this->segments_.perforations(),
+                                                 well_state_copy);
+
+        well_copy.calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
+        const double dt = ebos_simulator.timeStepSize();
+        // solve equations
+        const bool converged = well_copy.iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
+                                                                    deferred_logger);
+
+        // fetch potentials (sign is updated on the outside).
+        well_potentials.clear();
+        well_potentials.resize(np, 0.0);
+        for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
+            const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
+            well_potentials[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
+        }
+        debug_cost_counter_ += well_copy.debug_cost_counter_;
+        return converged;
+    }
 
     template <typename TypeTag>
     void
@@ -1227,6 +1300,13 @@ namespace Opm
         }
     }
 
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    updateIPRImplicit(const Simulator& ebos_simulator, WellState& well_state, DeferredLogger& deferred_logger)
+    {
+    }    
+
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
@@ -1441,7 +1521,8 @@ namespace Opm
                              const Well::ProductionControls& prod_controls,
                              WellState& well_state,
                              const GroupState& group_state,
-                             DeferredLogger& deferred_logger)
+                             DeferredLogger& deferred_logger, 
+                             const bool allow_switch /*true*/)
     {
         //if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
 
@@ -1586,9 +1667,9 @@ namespace Opm
                 } else {
                     this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
                 }
-                // We reset the well status to it's original state. Status is updated
+                // We reset the well status to it's original state. Status is updated 
                 // on the outside based on operability status
-                this->wellStatus_ = well_status;
+                // this->wellStatus_ = well_status;
             }
             std::string message = fmt::format("   Well {} converged in {} inner iterations ("
                                                     "{} control/status switches).", this->name(), it, switch_count);

opm/simulators/wells/SingleWellState.cpp

@@ -42,6 +42,8 @@ SingleWellState::SingleWellState(const std::string& name_,
     , temperature(temp)
     , well_potentials(pu_.num_phases)
     , productivity_index(pu_.num_phases)
+    , ipr_a(pu_.num_phases)
+    , ipr_b(pu_.num_phases)
     , surface_rates(pu_.num_phases)
     , reservoir_rates(pu_.num_phases)
     , prev_surface_rates(pu_.num_phases)
@@ -89,6 +91,8 @@ void SingleWellState::shut() {
     std::fill(this->prev_surface_rates.begin(), this->prev_surface_rates.end(), 0);
     std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0);
     std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
+    std::fill(this->ipr_a.begin(), this->ipr_a.end(), 0);
+    std::fill(this->ipr_b.begin(), this->ipr_b.end(), 0);
 
     auto& connpi = this->perf_data.prod_index;
     connpi.assign(connpi.size(), 0);
@@ -305,6 +309,8 @@ bool SingleWellState::operator==(const SingleWellState& rhs) const
            this->phase_mixing_rates == rhs.phase_mixing_rates &&
            this->well_potentials == rhs.well_potentials &&
            this->productivity_index == rhs.productivity_index &&
+           this->ipr_a == rhs.ipr_a &&
+           this->ipr_b == rhs.ipr_b &&
            this->surface_rates == rhs.surface_rates &&
            this->reservoir_rates == rhs.reservoir_rates &&
            this->prev_surface_rates == rhs.prev_surface_rates &&

opm/simulators/wells/SingleWellState.hpp

@@ -61,6 +61,8 @@ public:
         serializer(phase_mixing_rates);
         serializer(well_potentials);
         serializer(productivity_index);
+        serializer(ipr_a);
+        serializer(ipr_b);        
         serializer(surface_rates);
         serializer(reservoir_rates);
         serializer(prev_surface_rates);
@@ -98,6 +100,8 @@ public:
 
     std::vector<double> well_potentials;
     std::vector<double> productivity_index;
+    std::vector<double> ipr_a;
+    std::vector<double> ipr_b;    
     std::vector<double> surface_rates;
     std::vector<double> reservoir_rates;
     std::vector<double> prev_surface_rates;

opm/simulators/wells/StandardWell.hpp

@@ -212,7 +212,8 @@ namespace Opm
                                         const Well::ProductionControls& prod_controls,
                                         WellState& well_state,
                                         const GroupState& group_state,
-                                        DeferredLogger& deferred_logger) override;
+                                        DeferredLogger& deferred_logger, 
+                                        const bool allow_switch = true) override;
 
         /// \brief Wether the Jacobian will also have well contributions in it.
         virtual bool jacobianContainsWellContributions() const override
@@ -240,6 +241,8 @@ namespace Opm
             const double alq_value,
             DeferredLogger& deferred_logger) const override;
 
+        void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger);                
+
         virtual void computeWellRatesWithBhp(
             const Simulator& ebosSimulator,
             const double& bhp,
@@ -321,6 +324,10 @@ namespace Opm
             DeferredLogger& deferred_logger,
             const WellState &well_state) const;
 
+        bool computeWellPotentialsImplicit(const Simulator& ebos_simulator,
+                                           std::vector<double>& well_potentials,
+                                           DeferredLogger& deferred_logger) const;               
+
         virtual double getRefDensity() const override;
 
         // get the mobility for specific perforation

opm/simulators/wells/StandardWell_impl.hpp

@@ -832,6 +832,89 @@ namespace Opm
         this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
     }
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateIPRImplicit(const Simulator& ebosSimulator,
+                      WellState& well_state, 
+                      DeferredLogger& deferred_logger)
+    {   
+        // Compute IPR based on *converged* well-equation:
+        // For a component rate r the derivative dr/dbhp is obtained by 
+        // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial control_value)
+        // where Eq(x)=0 is the well equation setup with bhp control and primary varables x 
+        //StandardWell<TypeTag> well_copy(*this);
+
+        // We shouldn't have zero rates at this stage, but check
+        bool zero_rates;
+        auto rates = well_state.well(this->index_of_well_).surface_rates;
+        zero_rates = true;
+        for (std::size_t p = 0; p < rates.size(); ++p) {
+            zero_rates &= rates[p] == 0.0;
+        }
+        auto& ws = well_state.well(this->index_of_well_);
+        if (zero_rates) {
+            const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, reverting to explicit IPR-calulations", this->name());
+            deferred_logger.debug(msg);
+            updateIPR(ebosSimulator, deferred_logger);
+            for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
+                const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
+                ws.ipr_a[idx] = this->ipr_a_[comp_idx];
+                ws.ipr_b[idx] = this->ipr_b_[comp_idx];
+            }
+        } else {
+            const auto& group_state  = ebosSimulator.problem().wellModel().groupState();
+
+            // XXX maybe don't update this
+            std::fill(ws.ipr_a.begin(), ws.ipr_a.end(), 0.);
+            std::fill(ws.ipr_b.begin(), ws.ipr_b.end(), 0.);
+            //WellState well_state_copy = well_state;    
+            auto inj_controls = Well::InjectionControls(0);
+            auto prod_controls = Well::ProductionControls(0);
+            prod_controls.addControl(Well::ProducerCMode::BHP);
+            prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
+
+            //  Set current control to bhp, and bhp value in state, modify bhp limit in control object.
+            const auto cmode = ws.production_cmode;
+            ws.production_cmode = Well::ProducerCMode::BHP;
+            const double dt = ebosSimulator.timeStepSize();
+            assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+
+            const double nEq = this->primary_variables_.numWellEq();
+            BVectorWell rhs(1);
+            rhs[0].resize(nEq);
+            // rhs = 0 except -1 for control eq
+            for (size_t i=0; i < nEq; ++i){
+                rhs[0][i] = 0.0;            
+            }
+            rhs[0][Bhp] = -1.0;
+
+            BVectorWell x_well(1);
+            x_well[0].resize(nEq);
+            this->linSys_.solve(rhs, x_well);
+
+            for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
+                EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
+                const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
+                for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
+                    ws.ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
+                }
+                ws.ipr_a[idx] = ws.ipr_b[idx]*ws.bhp - comp_rate.value();
+                
+                //for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
+                //    this->ipr_b_[comp_idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
+                //}
+                // XXX maybe don't update this 
+                //this->ipr_a_[comp_idx] = this->ipr_b_[comp_idx]*ws.bhp - comp_rate.value();
+                // For ipr in well_state use same ordering as potentials etc.
+                //const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
+                //ws.ipr_a[idx] = this->ipr_a_[comp_idx];
+                //ws.ipr_b[idx] = this->ipr_b_[comp_idx];
+            }
+            // reset cmode
+            ws.production_cmode = cmode;
+        }
+    }
 
     template<typename TypeTag>
     void
@@ -1492,6 +1575,63 @@ namespace Opm
         return potentials;
     }
 
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    computeWellPotentialsImplicit(const Simulator& ebos_simulator,
+                                  std::vector<double>& well_potentials,
+                                  DeferredLogger& deferred_logger) const
+    {
+        // Create a copy of the well.
+        // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials 
+        // is allready a copy, but not from other calls.
+        StandardWell<TypeTag> well_copy(*this);
+
+        // store a copy of the well state, we don't want to update the real well state
+        WellState well_state_copy = ebos_simulator.problem().wellModel().wellState();
+        const auto& group_state = ebos_simulator.problem().wellModel().groupState();
+        auto& ws = well_state_copy.well(this->index_of_well_);
+
+        // get current controls
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        auto inj_controls = well_copy.well_ecl_.isInjector()
+            ? well_copy.well_ecl_.injectionControls(summary_state)
+            : Well::InjectionControls(0);
+        auto prod_controls = well_copy.well_ecl_.isProducer()
+            ? well_copy.well_ecl_.productionControls(summary_state) :
+            Well::ProductionControls(0);
+
+        // prepare/modify well state and control
+        well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
+        
+        // initialize rates from previous potentials
+        const int np = this->number_of_phases_;
+        bool trivial = true;
+        for (int phase = 0; phase < np; ++phase){
+            trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
+        }
+        if (!trivial) {
+            const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
+            for (int phase = 0; phase < np; ++phase) {
+                ws.surface_rates[phase] = sign * ws.well_potentials[phase];
+            }
+        }
+
+        well_copy.calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
+        const double dt = ebos_simulator.timeStepSize();
+        // iterate to get a solution at the given bhp.
+        const bool converged = well_copy.iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
+                                           deferred_logger);
+
+        // fetch potentials (sign is updated on the outside).
+        well_potentials.clear();
+        well_potentials.resize(np, 0.0);
+        for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
+            const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
+            well_potentials[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
+        }
+        return converged;
+    }
 
 
     template<typename TypeTag>
@@ -1554,29 +1694,35 @@ namespace Opm
             return;
         }
 
-        // does the well have a THP related constraint?
+        bool converged_implicit = false;
-        const auto& summaryState = ebosSimulator.vanguard().summaryState();
+        if (this->param_.local_well_solver_control_switching_) {
-        if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
+            converged_implicit = computeWellPotentialsImplicit(ebosSimulator, well_potentials, deferred_logger);
-            // get the bhp value based on the bhp constraints
+        }
-            double bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
+        if (!converged_implicit) {        
+            // does the well have a THP related constraint?
+            const auto& summaryState = ebosSimulator.vanguard().summaryState();
+            if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
+                // get the bhp value based on the bhp constraints
+                double bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
 
-            // In some very special cases the bhp pressure target are
+                // In some very special cases the bhp pressure target are
-            // temporary violated. This may lead to too small or negative potentials
+                // temporary violated. This may lead to too small or negative potentials
-            // that could lead to premature shutting of wells.
+                // that could lead to premature shutting of wells.
-            // As a remedy the bhp that gives the largest potential is used.
+                // As a remedy the bhp that gives the largest potential is used.
-            // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp,
+                // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp,
-            // and the potentials will be computed using the limit as expected.
+                // and the potentials will be computed using the limit as expected.
-            const auto& ws = well_state.well(this->index_of_well_);
+                const auto& ws = well_state.well(this->index_of_well_);
-            if (this->isInjector())
+                if (this->isInjector())
-                bhp = std::max(ws.bhp, bhp);
+                    bhp = std::max(ws.bhp, bhp);
-            else
+                else
-                bhp = std::min(ws.bhp, bhp);
+                    bhp = std::min(ws.bhp, bhp);
 
-            assert(std::abs(bhp) != std::numeric_limits<double>::max());
+                assert(std::abs(bhp) != std::numeric_limits<double>::max());
-            computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
+                computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
-        } else {
+            } else {
-            // the well has a THP related constraint
+                // the well has a THP related constraint
-            well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
+                well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
+            }
         }
 
         this->checkNegativeWellPotentials(well_potentials,
@@ -2174,32 +2320,40 @@ namespace Opm
                                const Well::ProductionControls& prod_controls,
                                WellState& well_state,
                                const GroupState& group_state,
-                               DeferredLogger& deferred_logger)
+                               DeferredLogger& deferred_logger, 
+                               const bool allow_switch /*true*/)
     {
         const int max_iter = this->param_.max_inner_iter_wells_;
-
+        
         int it = 0;
         bool converged;
         bool relax_convergence = false;
         this->regularize_ = false;
         const auto& summary_state = ebosSimulator.vanguard().summaryState();
 
-        // Max status switch frequency should be 2 to avoid getting stuck in cycle
+        // Max status switch frequency should be 2 to avoid getting stuck in cycle 
-        constexpr int min_its_after_switch = 2;
+        constexpr int min_its_after_switch = 4;
         int its_since_last_switch = min_its_after_switch;
         int switch_count= 0;
-        const auto well_status = this->wellStatus_;
+        // if we fail to solve eqs, we reset status/control before leaving
+        const auto well_status_orig = this->wellStatus_;
+        auto well_status_cur = well_status_orig;
+        int status_switch_count = 0;
         const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
         bool changed = false;
-        bool final_check = false;
+        bool final_check = false; 
         do {
             its_since_last_switch++;
             if (allow_switching && its_since_last_switch >= min_its_after_switch){
                 const double wqTotal = this->primary_variables_.eval(WQTotal).value();
-                changed = this->updateWellControlAndStatusLocalIteration(ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger);
+                changed = this->updateWellControlAndStatusLocalIteration(ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger, allow_switch); 
                 if (changed){
                     its_since_last_switch = 0;
                     switch_count++;
+                    if (well_status_cur != this->wellStatus_) {
+                        well_status_cur = this->wellStatus_;
+                        status_switch_count++;
+                    }
                 }
                 if (!changed && final_check) {
                     break;
@@ -2207,7 +2361,7 @@ namespace Opm
                     final_check = false;
                 }
             }
-
+  
             assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
 
             if (it > this->param_.strict_inner_iter_wells_) {
@@ -2226,14 +2380,20 @@ namespace Opm
                     its_since_last_switch = min_its_after_switch;
                 } else {
                     break;
-                }
+                }   
             }
 
             ++it;
             solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
             initPrimaryVariablesEvaluation();
-        } while (it < max_iter);
 
+            if (status_switch_count > 10) {
+                // constantly changing status, give up (do operability check on the outside)
+                converged = false;
+                break;
+            }
+        } while (it < max_iter);
+        
         if (converged) {
             if (allow_switching){
                 // update operability if status change
@@ -2251,11 +2411,12 @@ namespace Opm
                 // For this to happen, isOperableAndSolvable() must change from true to false,
                 // and (until the most recent commit) the well needs to be open for this to trigger.
                 // Hence, the resetting of status.
-                this->wellStatus_ = well_status;
+                //this->wellStatus_ = well_status;
             }
         } else {
+            this->wellStatus_ = well_status_orig;
             const std::string message = fmt::format("   Well {} did not converged in {} inner iterations ("
-                                                    "{} control/status switches).", this->name(), it, switch_count);
+                                                    "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
             deferred_logger.debug(message);
             // add operability here as well ?
         }

opm/simulators/wells/VFPHelpers.cpp

@@ -499,6 +499,90 @@ double findTHP(const std::vector<double>& bhp_array,
     return thp;
 }
 
+std::pair<double, double> 
+getMinimumBHPCoordinate(const VFPProdTable& table,
+                        const double thp,
+                        const double wfr,
+                        const double gfr,
+                        const double alq)
+{   
+    double flo_at_bhp_min = 0.0; // start by checking flo=0
+    auto flo_i = detail::findInterpData(flo_at_bhp_min, table.getFloAxis());
+    auto thp_i = detail::findInterpData( thp, table.getTHPAxis());
+    auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis());
+    auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis());
+    auto alq_i = detail::findInterpData( alq, table.getALQAxis());
+
+    detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+    double bhp_min = bhp_i.value;
+    std::vector<double> flos = table.getFloAxis();
+    for (size_t i = 0; i < flos.size(); ++i) {
+        flo_i = detail::findInterpData(flos[i], table.getFloAxis());
+        bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+        if (bhp_i.value < bhp_min){
+            bhp_min = bhp_i.value;
+            flo_at_bhp_min = flos[i];
+        }
+    }
+    // return negative flo
+    return std::make_pair(-flo_at_bhp_min, bhp_min); 
+}      
+
+std::optional<std::pair<double, double>> 
+intersectWithIPR(const VFPProdTable& table,
+                 const double thp,
+                 const double wfr,
+                 const double gfr,
+                 const double alq, 
+                 const double ipr_a,
+                 const double ipr_b)
+{   
+    // NOTE: ipr-line is q=b*bhp - a!
+    // ipr is given for negative flo, so
+    // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b  
+    auto thp_i = detail::findInterpData( thp, table.getTHPAxis());
+    auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis());
+    auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis());
+    auto alq_i = detail::findInterpData( alq, table.getALQAxis());
+
+    if (ipr_b == 0.0) {
+        // this shouldn't happen, but deal with it to be safe
+        auto flo_i = detail::findInterpData(-ipr_a, table.getFloAxis());
+        detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+        return std::make_pair(ipr_b, bhp_i.value);
+    }
+    // find largest flo (flo_x) for which y = bhp(flo) + (flo-a)/b = 0 and dy/dflo > 0
+    double flo_x = -1.0;
+    double flo0, flo1;
+    double y0, y1;
+    flo0 = 0.0; // start by checking flo=0
+    auto flo_i = detail::findInterpData(flo0, table.getFloAxis());
+    detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+    y0 = bhp_i.value - ipr_a/ipr_b; // +0.0/ipr_b
+
+    std::vector<double> flos = table.getFloAxis();
+    for (size_t i = 0; i < flos.size(); ++i) {
+        flo1 = flos[i];
+        flo_i = detail::findInterpData(flo1, table.getFloAxis());
+        bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+        y1 = bhp_i.value + (flo1 - ipr_a)/ipr_b;
+        if (y0 < 0 && y1 >= 0){
+            // crossing with positive slope
+            double w = -y0/(y1-y0);
+            w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1)
+            flo_x = flo0 + w*(flo1 - flo0);
+        }
+        flo0 = flo1;
+        y0 = y1;
+    }
+    // return (last) intersection if found (negative flo)
+    if (flo_x >= 0) {
+        return std::make_pair(-flo_x, -(flo_x - ipr_a)/ipr_b);
+    } else {
+        return std::nullopt;
+    }
+} 
+
 template <typename T>
 T getFlo(const VFPProdTable& table,
          const T& aqua,

opm/simulators/wells/VFPHelpers.hpp

@@ -25,6 +25,7 @@
 #include <functional>
 #include <map>
 #include <vector>
+#include <optional>
 
 /**
  * This file contains a set of helper functions used by VFPProd / VFPInj.
@@ -181,6 +182,28 @@ double findTHP(const std::vector<double>& bhp_array,
                const std::vector<double>& thp_array,
                double bhp);
 
+/**
+* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq
+*/
+std::pair<double, double> 
+getMinimumBHPCoordinate(const VFPProdTable& table,
+                        const double thp,
+                        const double wfr,
+                        const double gfr,
+                        const double alq);
+
+/**
+* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection
+* if it exists
+*/  
+std::optional<std::pair<double, double>> 
+intersectWithIPR(const VFPProdTable& table,
+                 const double thp,
+                 const double wfr,
+                 const double gfr,
+                 const double alq, 
+                 const double ipr_a,
+                 const double ipr_b);                        
 
 } // namespace detail
 

opm/simulators/wells/VFPProdProperties.cpp

@@ -137,6 +137,22 @@ bhpwithflo(const std::vector<double>& flos,
     return bhps;
 }
 
+double
+VFPProdProperties::
+minimumBHP(const int table_id,
+           const double thp,
+           const double wfr,
+           const double gfr,
+           const double alq) const
+{
+    // Get the table
+    const VFPProdTable& table = detail::getTable(m_tables, table_id);
+    const auto retval = detail::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq);
+    // returned pair is (flo, bhp)
+    return retval.second;
+}
+
+
 
 void VFPProdProperties::addTable(const VFPProdTable& new_table) {
     this->m_tables.emplace( new_table.getTableNum(), new_table );
@@ -176,7 +192,8 @@ EvalWell VFPProdProperties::bhp(const int table_id,
 
     detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
 
-    bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
+   //bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
+   bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo* flo);
 
     bhp.setValue(bhp_val.value);
     return bhp;

opm/simulators/wells/VFPProdProperties.hpp

@@ -129,7 +129,11 @@ public:
         return m_tables.empty();
     }
 
-
+    /**
+     * Returns minimum bhp for given thp, wfr, gfr and alq 
+     */
+    double minimumBHP(const int table_id, const double thp, 
+                      const double wfr, const double gfr, const double alq) const;
 protected:
     // calculate a group bhp values with a group of flo rate values
     std::vector<double> bhpwithflo(const std::vector<double>& flos,

opm/simulators/wells/VFPProperties.hpp

@@ -92,6 +92,23 @@ public:
         return detail::getGFR(table, aqua, liquid, vapour);
     }
 
+    std::pair<double, double> 
+    getFloIPR(const int table_id, const std::size_t well_index) const {
+        std::pair<double,double>retval(0.0, 0.0);
+        const VFPProdTable& table = this->m_prod.getTable(table_id);
+        const auto& pu = well_state_.phaseUsage();
+        const auto& ipr_a= well_state_.well(well_index).ipr_a;
+        const auto& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
+        const auto& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;
+        const auto& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]]:0.0;
+        const auto& ipr_b= well_state_.well(well_index).ipr_b;
+        const auto& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
+        const auto& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;
+        const auto& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]]:0.0;
+        return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), 
+                                detail::getFlo(table, aqua_b, liquid_b, vapour_b));
+    }
+
 private:
     VFPInjProperties m_inj;
     VFPProdProperties m_prod;

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -387,6 +387,32 @@ calculateBhpFromThp(const WellState& well_state,
     return bhp_tab - dp_hydro + bhp_adjustment;
 }
 
+double
+WellBhpThpCalculator::
+calculateMinimumBhpFromThp(const WellState& well_state,
+                           const Well& well,
+                           const SummaryState& summaryState,
+                           const double rho) const
+{
+    assert(well_.isProducer()); // only producers can go here for now
+
+    const double thp_limit = well_.getTHPConstraint(summaryState);
+    
+    const auto& controls = well.productionControls(summaryState);
+    const auto& wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
+    const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
+
+    const double bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number,
+                                                                        thp_limit, wfr, gfr, 
+                                                                        well_.getALQ(well_state));
+ 
+    const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
+    const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit);
+    const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, 
+                                                                      rho, well_.gravity());
+    return bhp_min - dp_hydro + bhp_adjustment;
+}
+
 double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const
 {
     return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit);
@@ -839,6 +865,96 @@ bruteForceBracket(const std::function<double(const double)>& eq,
     return bracket_found;
 }
 
+bool WellBhpThpCalculator::
+isStableSolution(const WellState& well_state,
+                 const Well& well,
+                 const std::vector<double>& rates,
+                 const SummaryState& summaryState,
+                 DeferredLogger& deferred_logger) const
+{
+    assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
+    assert(well_.isProducer()); // only valid for producers 
+
+    static constexpr int Gas = BlackoilPhases::Vapour;
+    static constexpr int Oil = BlackoilPhases::Liquid;
+    static constexpr int Water = BlackoilPhases::Aqua;
+
+    const double aqua = rates[Water];
+    const double liquid = rates[Oil];
+    const double vapour = rates[Gas];
+    const double thp = well_.getTHPConstraint(summaryState);
+
+    const auto& controls = well.productionControls(summaryState);
+    const auto& wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
+    const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
+
+    const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
+    const bool use_vfpexplicit = well_.useVfpExplicit();
+
+    const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
+    const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp);
+    // XXX this needs to be fixed 
+    //assert(bhp_adjustment == 0.0);
+
+    const detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit);
+
+    if (bhp.dflo >= 0) {
+        return true;
+    } else {    // maybe check if ipr is available
+        const auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell());
+        return bhp.dflo + 1/ipr.second >= 0; 
+    }                  
+}
+
+std::optional<double> WellBhpThpCalculator::
+estimateStableBhp(const WellState& well_state,
+                  const Well& well,
+                  const std::vector<double>& rates,
+                  const double rho,
+                  const SummaryState& summaryState,
+                  DeferredLogger& deferred_logger) const
+{   
+    // Given a *converged* well_state with ipr, estimate bhp of the stable solution 
+    const auto& controls = well.productionControls(summaryState);
+    const double thp = well_.getTHPConstraint(summaryState);
+    const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
+
+    const double aqua = rates[BlackoilPhases::Aqua];
+    const double liquid = rates[BlackoilPhases::Liquid];
+    const double vapour = rates[BlackoilPhases::Vapour];
+    double flo = detail::getFlo(table, aqua, liquid, vapour);
+    double wfr, gfr;
+    if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) {
+        wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
+        gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
+    } else {
+        wfr = detail::getWFR(table, aqua, liquid, vapour);  
+        gfr = detail::getGFR(table, aqua, liquid, vapour);   
+    }
+    if (wfr <= 0.0 && gfr <= 0.0) {
+        // warning message
+    }
+    auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell());
+    
+    const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
+    const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp);
+    // XXX this needs to be fixed 
+    //assert(bhp_adjustment == 0.0);
+
+    const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, 
+                                                                      rho, well_.gravity());
+    if (ipr.first <= 0.0) {
+        // error message
+    }
+    const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first+ipr.second*dp_hydro, ipr.second);
+    if (retval.has_value()) {
+        // returned pair is (flo, bhp)
+        return retval.value().second - dp_hydro;
+    } else {
+        return std::nullopt;
+    }
+}        
+
 #define INSTANCE(...) \
 template __VA_ARGS__ WellBhpThpCalculator:: \
 calculateBhpFromThp<__VA_ARGS__>(const WellState&, \

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -98,6 +98,25 @@ public:
                                const double rho,
                                DeferredLogger& deferred_logger) const;
 
+  double calculateMinimumBhpFromThp(const WellState& well_state,
+                                    const Well& well,
+                                    const SummaryState& summaryState,
+                                    const double rho) const;      
+
+  bool isStableSolution(const WellState& well_state,
+                        const Well& well,
+                        const std::vector<double>& rates,
+                        const SummaryState& summaryState,
+                        DeferredLogger& deferred_logger) const;   
+
+  std::optional<double>
+  estimateStableBhp (const WellState& well_state,
+                    const Well& well,
+                    const std::vector<double>& rates,
+                    const double rho,
+                    const SummaryState& summaryState,
+                    DeferredLogger& deferred_logger) const;                                                                             
+
 private:
     //! \brief Compute BHP from THP limit for an injector - implementation.
     template<class ErrorPolicy>

opm/simulators/wells/WellInterface.hpp

@@ -247,7 +247,8 @@ public:
                                                   const Well::InjectionControls& inj_controls,
                                                   const Well::ProductionControls& prod_controls,
                                                   const double WQTotal,
-                                                  DeferredLogger& deferred_logger);
+                                                  DeferredLogger& deferred_logger, 
+                                                  const bool allow_switch=true);
 
     virtual void updatePrimaryVariables(const SummaryState& summary_state,
                                         const WellState& well_state,
@@ -414,7 +415,12 @@ protected:
                                             const WellProductionControls& prod_controls,
                                             WellState& well_state,
                                             const GroupState& group_state,
-                                            DeferredLogger& deferred_logger) = 0;
+                                            DeferredLogger& deferred_logger, 
+                                            const bool allow_switch = true) = 0;
+
+    virtual void updateIPRImplicit(const Simulator& ebosSimulator,
+                                   WellState& well_state,
+                                   DeferredLogger& deferred_logger) = 0;                                            
 
     bool iterateWellEquations(const Simulator& ebosSimulator,
                               const double dt,
@@ -422,6 +428,27 @@ protected:
                               const GroupState& group_state,
                               DeferredLogger& deferred_logger);
 
+    bool robustWellSolveWithTHP(const Simulator& ebos_simulator,
+                                const double dt,
+                                const Well::InjectionControls& inj_controls,
+                                const Well::ProductionControls& prod_controls,                                
+                                WellState& well_state,
+                                const GroupState& group_state,
+                                DeferredLogger& deferred_logger);
+
+    std::optional<double> estimateOperableBhp(const Simulator& ebos_simulator,
+                                              const double dt,
+                                              WellState& well_state,
+                                              const GroupState& group_state,
+                                              const SummaryState& summary_state,
+                                              DeferredLogger& deferred_logger);        
+
+    bool solveWellWithBhp(const Simulator& ebos_simulator,
+                          const double dt,
+                          const double bhp,
+                          WellState& well_state,
+                          DeferredLogger& deferred_logger);                                                                                       
+
     bool solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
                              DeferredLogger& deferred_logger);
 

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -744,4 +744,47 @@ checkNegativeWellPotentials(std::vector<double>& well_potentials,
     }
 }
 
+void WellInterfaceGeneric::
+prepareForPotentialCalculations(const SummaryState& summary_state,
+                                WellState& well_state, 
+                                Well::InjectionControls& inj_controls,
+                                Well::ProductionControls& prod_controls)
+{
+    const bool has_thp = this->wellHasTHPConstraints(summary_state);
+    auto& ws = well_state.well(this->index_of_well_);
+    // Modify control (only pressure constraints) and set new target if needed.
+    // Also set value for current target in state
+    if (this->isInjector()) {
+        inj_controls.clearControls();
+        inj_controls.addControl(Well::InjectorCMode::BHP);
+        if (has_thp){
+            inj_controls.addControl(Well::InjectorCMode::THP);
+        }
+        if (!(ws.injection_cmode == Well::InjectorCMode::BHP)){
+            if (has_thp){
+                ws.injection_cmode = Well::InjectorCMode::THP;
+                ws.thp = this->getTHPConstraint(summary_state);
+            } else {
+                ws.injection_cmode = Well::InjectorCMode::BHP;
+                ws.bhp = inj_controls.bhp_limit;
+            }
+        } 
+    } else { // producer
+        prod_controls.clearControls();
+        prod_controls.addControl(Well::ProducerCMode::BHP);
+        if (has_thp){
+            prod_controls.addControl(Well::ProducerCMode::THP);
+        }
+        if (!(ws.production_cmode == Well::ProducerCMode::BHP)){
+            if (has_thp){
+                ws.production_cmode = Well::ProducerCMode::THP;
+                ws.thp = this->getTHPConstraint(summary_state);
+            } else {
+                ws.production_cmode = Well::ProducerCMode::BHP;
+                ws.bhp = prod_controls.bhp_limit;
+            }
+        } 
+    }    
+}
+
 } // namespace Opm

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -245,6 +245,11 @@ protected:
                                      const bool checkOperability,
                                      DeferredLogger& deferred_logger);
 
+    void prepareForPotentialCalculations(const SummaryState& summary_state,
+                                         WellState& well_state, 
+                                         Well::InjectionControls& inj_controls,
+                                         Well::ProductionControls& prod_controls);                                  
+
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {
         bool isOperableAndSolvable() const {

opm/simulators/wells/WellInterface_impl.hpp

@@ -264,11 +264,12 @@ namespace Opm
                                              const Well::InjectionControls& inj_controls,
                                              const Well::ProductionControls& prod_controls,
                                              const double wqTotal,
-                                             DeferredLogger& deferred_logger)
+                                             DeferredLogger& deferred_logger, 
+                                             const bool allow_switch /*true*/)
     {
         const auto& summary_state = ebos_simulator.vanguard().summaryState();
         const auto& schedule = ebos_simulator.vanguard().schedule();
-
+        
         if (this->wellUnderZeroRateTarget(summary_state, well_state) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
            return false;
         }
@@ -276,30 +277,33 @@ namespace Opm
         const double sgn = this->isInjector() ? 1.0 : -1.0;
         if (!this->wellIsStopped()){
             if (wqTotal*sgn <= 0.0){
+                //std::cout << "Stopping well:" << this->name() << std::endl;
                 this->stopWell();
                 return true;
             } else {
                 bool changed = false;
-                auto& ws = well_state.well(this->index_of_well_);
+                if (allow_switch) {
-                const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
+                    auto& ws = well_state.well(this->index_of_well_);
-                                                                  prod_controls.hasControl(Well::ProducerCMode::GRUP);
+                    const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
+                                                                    prod_controls.hasControl(Well::ProducerCMode::GRUP);
 
-                changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
+                    changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
-                // TODO: with current way, the checkGroupConstraints might overwrite the result from checkIndividualConstraints, which remains to be investigated
+                    // TODO: with current way, the checkGroupConstraints might overwrite the result from checkIndividualConstraints, which remains to be investigated
-                if (hasGroupControl) {
+                    if (hasGroupControl) {
-                    changed = this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
+                        changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
-                }
+                    }
-
+
-                if (changed) {
+                    if (changed) {
-                    const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
+                        const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
-                                                                     ws.production_cmode == Well::ProducerCMode::THP;
+                                                                        ws.production_cmode == Well::ProducerCMode::THP;
-                    if (!thp_controlled){
+                        if (!thp_controlled){
-                        // don't call for thp since this might trigger additional local solve
+                            // don't call for thp since this might trigger additional local solve
-                        updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
+                            updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
-                    } else {
+                        } else {
-                        ws.thp = this->getTHPConstraint(summary_state);
+                            ws.thp = this->getTHPConstraint(summary_state);
+                        }  
+                        updatePrimaryVariables(summary_state, well_state, deferred_logger);
                     }
-                    updatePrimaryVariables(summary_state, well_state, deferred_logger);
                 }
                 return changed;
             }
@@ -311,20 +315,29 @@ namespace Opm
             const bool has_thp = this->wellHasTHPConstraints(summary_state);
             if (has_thp){
                 // calculate bhp from thp-limit (using explicit fractions zince zero rate)
-                // TODO: this will often be too strict condition for re-opening, a better
+                // TODO: this will often be too strict condition for re-opening, a better 
-                // option is probably minimum bhp on current vfp-curve, but some more functionality
+                // option is probably minimum bhp on current vfp-curve, but some more functionality 
-                // is needed for this option to be robustly implemented.
+                // is needed for this option to be robustly implemented. 
                 std::vector<double> rates(this->num_components_);
-                const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
+                //const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
                 if (this->isInjector()){
+                    const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
                     inj_limit = std::min(bhp_thp, inj_controls.bhp_limit);
                 } else {
-                    prod_limit = std::max(bhp_thp, prod_controls.bhp_limit);
+                    const double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
+                    prod_limit = std::max(bhp_min, prod_controls.bhp_limit);
+                    //auto prates = well_state.well(this->index_of_well_).prev_surface_rates;
+                    //std::cout << this->name() << ": Min bhp: " << bhp_min << "  prod limit: " << prod_limit << " prev rates: " << prates[0] << " " << prates[1] << " " << prates[2] << std::endl; 
                 }
             }
             const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
             if (bhp_diff > 0){
+                //std::cout << "Re-opening well:" << this->name() << std::endl; 
                 this->openWell();
+                well_state.well(this->index_of_well_).bhp = prod_limit;
+                if (has_thp) {
+                    well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
+                }
                 return true;
             } else {
                 return false;
@@ -346,7 +359,9 @@ namespace Opm
 
         WellState well_state_copy = well_state;
         auto& ws = well_state_copy.well(this->indexOfWell());
-
+        if (ws.production_cmode == Well::ProducerCMode::GRUP) {
+            ws.production_cmode = Well::ProducerCMode::BHP;
+        }
         updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
         calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
         const auto& summary_state = simulator.vanguard().summaryState();
@@ -354,12 +369,7 @@ namespace Opm
         initPrimaryVariablesEvaluation();
 
         if (this->isProducer()) {
-            const auto& schedule = simulator.vanguard().schedule();
+            gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
-            const auto report_step = simulator.episodeIndex();
-            const auto& glo = schedule.glo(report_step);
-            if (glo.active()) {
-                gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
-            }
         }
 
         WellTestState welltest_state_temp;
@@ -448,9 +458,13 @@ namespace Opm
             if (!this->param_.local_well_solver_control_switching_){
                 converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
             } else {
-                converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+                if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
+                    converged = robustWellSolveWithTHP(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+                } else {
+                    converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+                }
             }
-
+                 
         } catch (NumericalProblem& e ) {
             const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
             deferred_logger.warning("INNER_ITERATION_FAILED", msg);
@@ -459,6 +473,141 @@ namespace Opm
         return converged;
     }
 
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    robustWellSolveWithTHP(const Simulator& ebos_simulator,
+                           const double dt,
+                           const Well::InjectionControls& inj_controls,
+                           const Well::ProductionControls& prod_controls,                           
+                           WellState& well_state,
+                           const GroupState& group_state,
+                           DeferredLogger& deferred_logger)
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        bool is_operable = true;
+        bool converged;
+        converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+        if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP) {
+            auto rates = well_state.well(this->index_of_well_).surface_rates;
+            this->adaptRatesForVFP(rates);
+            // FIX !!!!!!!!!!
+            // first, check just dbhp/dflo
+
+            // update ipr for stability-check
+            this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger);
+            bool is_stable =  WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state, deferred_logger);
+            if (!is_stable) {
+                // solution converged to an unstable point! 
+                //this->operability_status_.use_vfpexplicit = true;
+                auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger);
+                if (!bhp_stable.has_value()) {
+                    // well is not operable (no intersection)
+                    is_operable = false;
+                } else {
+                    // solve equations for bhp_stable (in attempt to get closer to actual solution)
+                    const bool converged_bhp = solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger);
+                    if (!converged_bhp) {
+                        // 
+                    } else {
+                        // re-solve well equations
+                        converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+                    }
+                }
+            }            
+        } else if (!converged) {
+            // Well did not converge, switch to explicit fractions
+            this->operability_status_.use_vfpexplicit = true;
+            //auto well_state_copy = well_state;
+            // 
+            auto bhp_target = estimateOperableBhp(ebos_simulator, dt, well_state, group_state, summary_state, deferred_logger);
+            if (!bhp_target.has_value()) {
+                // well can't operate using explicit fractions
+                is_operable = false;
+            } else {      
+                // first solve well equations for bhp = bhp_target which should give solution close to actual
+                converged = solveWellWithBhp(ebos_simulator, dt, bhp_target.value(), well_state, deferred_logger);
+                if (!converged) {
+                    // debug
+                } else {
+                    // re-solve well equations
+                    converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+                }
+            }
+        }
+        if (!is_operable) {
+            this->stopWell();
+            if (this->wellHasTHPConstraints(summary_state)){
+                this->operability_status_.can_obtain_bhp_with_thp_limit = false;
+                this->operability_status_.obey_thp_limit_under_bhp_limit = false;
+            } else {
+                this->operability_status_.operable_under_only_bhp_limit = false;
+            }
+        }
+        return converged;
+    }
+
+    template<typename TypeTag>
+    std::optional<double>
+    WellInterface<TypeTag>::
+    estimateOperableBhp(const Simulator& ebos_simulator,
+                        const double dt,
+                        WellState& well_state,
+                        const GroupState& group_state,
+                        const SummaryState& summary_state,
+                        DeferredLogger& deferred_logger)
+    {   
+        // Given an unconverged well or closed well, estimate an operable bhp (if any) 
+        // Get minimal bhp from vfp-curve
+        double bhp_min =  WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
+        // Solve
+        const bool converged = solveWellWithBhp(ebos_simulator, dt, bhp_min, well_state, deferred_logger);
+        if (!converged) {
+            //report problem - return null 
+        } else if (this->wellIsStopped()) {
+            return std::nullopt;
+        }
+        this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger);
+        auto rates = well_state.well(this->index_of_well_).surface_rates;
+        this->adaptRatesForVFP(rates);
+        return WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger);
+    } 
+
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    solveWellWithBhp(const Simulator& ebos_simulator,
+                     const double dt,
+                     const double bhp,
+                     WellState& well_state,
+                     DeferredLogger& deferred_logger)
+    {   
+        // Solve a well using single bhp-constraint (but close if not operable under this) 
+        auto group_state = GroupState(); // empty group
+        auto inj_controls = Well::InjectionControls(0);
+        auto prod_controls = Well::ProductionControls(0); 
+        auto& ws = well_state.well(this->index_of_well_);
+        auto cmode_inj = ws.injection_cmode;
+        auto cmode_prod = ws.production_cmode;
+        if (this->isInjector()) {
+            inj_controls.addControl(Well::InjectorCMode::BHP);
+            inj_controls.bhp_limit = bhp;
+            inj_controls.cmode = Well::InjectorCMode::BHP;
+            ws.injection_cmode = Well::InjectorCMode::BHP;  
+        } else {
+            prod_controls.addControl(Well::ProducerCMode::BHP);
+            prod_controls.bhp_limit = bhp;
+            prod_controls.cmode = Well::ProducerCMode::BHP;
+            ws.production_cmode = Well::ProducerCMode::BHP;  
+        }
+        // update well-state
+        ws.bhp = bhp;
+        // solve 
+        const bool converged =  this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*allow_switch*/false);
+        ws.injection_cmode = cmode_inj;
+        ws.production_cmode = cmode_prod;
+        return converged;
+    } 
 
     template<typename TypeTag>
     bool
@@ -717,8 +866,8 @@ namespace Opm
         const auto& schedule = ebos_simulator.vanguard().schedule();
         auto report_step_idx = ebos_simulator.episodeIndex();
         const auto& glo = schedule.glo(report_step_idx);
-        if(glo.active() && glo.has_well(well_name)) {
+        if(glo.has_well(well_name)) {
-            const auto increment = glo.gaslift_increment();
+            auto increment = glo.gaslift_increment();
             auto alq = well_state.getALQ(well_name);
             bool converged;
             while (alq > 0) {
@@ -751,8 +900,9 @@ namespace Opm
             deferred_logger.info(msg);
             return;
         }
+        const auto& well_ecl = this->wellEcl();
         const auto& schedule = ebos_simulator.vanguard().schedule();
-        const auto report_step_idx = ebos_simulator.episodeIndex();
+        auto report_step_idx = ebos_simulator.episodeIndex();
         const auto& glo = schedule.glo(report_step_idx);
         if (!glo.has_well(well_name)) {
             const std::string msg = fmt::format(
@@ -768,7 +918,6 @@ namespace Opm
             max_alq = *max_alq_optional;
         }
         else {
-            const auto& well_ecl = this->wellEcl();
             const auto& controls = well_ecl.productionControls(summary_state);
             const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
             const auto& alq_values = table.getALQAxis();
@@ -787,7 +936,7 @@ namespace Opm
     updateWellOperability(const Simulator& ebos_simulator,
                           const WellState& well_state,
                           DeferredLogger& deferred_logger)
-    {
+    {   
         if (this->param_.local_well_solver_control_switching_) {
             const bool success = updateWellOperabilityFromWellEq(ebos_simulator, well_state, deferred_logger);
             if (success) {
@@ -833,7 +982,7 @@ namespace Opm
         // equations should be converged at this stage, so only one it is needed
         bool converged = iterateWellEquations(ebos_simulator, dt, well_state_copy, group_state, deferred_logger);
         return converged;
-    }
+    }                          
 
     template<typename TypeTag>
     void