Merge pull request #4895 from GitPaean/steink-local_control_switch_option_0920_3

Option for doing well control/status switching during local well solve

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -174,6 +174,10 @@ template<class TypeTag, class MyTypeTag>
 struct UseAverageDensityMsWells {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct LocalWellSolveControlSwitching {
+    using type = UndefinedProperty;
+};
 // Network solver parameters
 template<class TypeTag, class MyTypeTag>
 struct NetworkMaxStrictIterations {
@@ -366,6 +370,11 @@ template<class TypeTag>
 struct UseAverageDensityMsWells<TypeTag, TTag::FlowModelParameters> {
     static constexpr bool value = false;
 };
+template<class TypeTag>
+struct LocalWellSolveControlSwitching<TypeTag, TTag::FlowModelParameters> {
+    static constexpr bool value = false;
+};
+
 // Network solver parameters
 template<class TypeTag>
 struct NetworkMaxStrictIterations<TypeTag, TTag::FlowModelParameters> {
@@ -530,6 +539,9 @@ namespace Opm
         /// Whether to approximate segment densities by averaging over segment and its outlet 
         bool use_average_density_ms_wells_;
 
+        /// Whether to allow control switching during local well solutions 
+        bool local_well_solver_control_switching_;
+
         /// Maximum number of iterations in the network solver before relaxing tolerance
         int network_max_strict_iterations_;
         
@@ -586,6 +598,7 @@ namespace Opm
             check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
             max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
             use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
+            local_well_solver_control_switching_ = EWOMS_GET_PARAM(TypeTag, bool, LocalWellSolveControlSwitching);
             nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver);
             std::string approach = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach);
             if (approach == "jacobi") {
@@ -651,6 +664,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet");
+            EWOMS_REGISTER_PARAM(TypeTag, bool, LocalWellSolveControlSwitching, "Allow control switching during local well solutions");
             EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxStrictIterations, "Maximum iterations in network solver before relaxing tolerance");
             EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxIterations, "Maximum number of iterations in the network solver before giving up");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, NonlinearSolver, "Choose nonlinear solver. Valid choices are newton or nldd.");

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -2258,6 +2258,9 @@ namespace Opm {
                     deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
                 }
             }
+            // If we're using local well solves that include control switches, they also update
+            // operability, so reset before main iterations begin
+            well->resetWellOperability();
         }
         updatePrimaryVariables(deferred_logger);
 

opm/simulators/wells/MultisegmentWell.hpp

@@ -256,6 +256,14 @@ namespace Opm
                                               const GroupState& group_state,
                                               DeferredLogger& deferred_logger) override;
 
+        virtual bool iterateWellEqWithSwitching(const Simulator& ebosSimulator,
+                                                const double dt,
+                                                const Well::InjectionControls& inj_controls,
+                                                const Well::ProductionControls& prod_controls,
+                                                WellState& well_state,
+                                                const GroupState& group_state,
+                                                DeferredLogger& deferred_logger) override;
+
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                                     const double dt,
                                                     const Well::InjectionControls& inj_controls,

opm/simulators/wells/MultisegmentWellGeneric.cpp

@@ -140,11 +140,9 @@ segmentNumberToIndex(const int segment_number) const
 template<typename Scalar>
 void
 MultisegmentWellGeneric<Scalar>::
-detectOscillations(const std::vector<double>& measure_history,
+detectOscillations(const std::vector<double>& measure_history, bool& oscillate, bool& stagnate) const
-                   const int it,
-                   bool& oscillate,
-                   bool& stagnate) const
 {
+    const auto it = measure_history.size() - 1;
     if ( it < 2 ) {
         oscillate = false;
         stagnate = false;

opm/simulators/wells/MultisegmentWellGeneric.hpp

@@ -65,7 +65,6 @@ protected:
 
     /// Detect oscillation or stagnation based on the residual measure history
     void detectOscillations(const std::vector<double>& measure_history,
-                            const int it,
                             bool& oscillate,
                             bool& stagnate) const;
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1339,7 +1339,7 @@ namespace Opm
             bool is_oscillate = false;
             bool is_stagnate = false;
 
-            this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
+            this->detectOscillations(measure_history, is_oscillate, is_stagnate);
             // TODO: maybe we should have more sophisticated strategy to recover the relaxation factor,
             // for example, to recover it to be bigger
 
@@ -1412,7 +1412,179 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    iterateWellEqWithSwitching(const Simulator& ebosSimulator,
+                             const double dt,
+                             const Well::InjectionControls& inj_controls,
+                             const Well::ProductionControls& prod_controls,
+                             WellState& well_state,
+                             const GroupState& group_state,
+                             DeferredLogger& deferred_logger)
+    {
+        //if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
 
+        const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
+
+        {
+            // getWellFiniteResiduals returns false for nan/inf residuals
+            const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+            if(!isFinite)
+                return false;
+        }
+
+        std::vector<std::vector<Scalar> > residual_history;
+        std::vector<double> measure_history;
+        int it = 0;
+        // relaxation factor
+        double relaxation_factor = 1.;
+        const double min_relaxation_factor = 0.6;
+        bool converged = false;
+        int stagnate_count = 0;
+        bool relax_convergence = false;
+        this->regularize_ = false;
+
+        // Max status switch frequency should be 2 to avoid getting stuck in cycle 
+        const int min_its_after_switch = 2;
+        int its_since_last_switch = min_its_after_switch;
+        int switch_count= 0;
+        const auto well_status = this->wellStatus_;
+        const auto& summary_state = ebosSimulator.vanguard().summaryState();
+        const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
+        bool changed = false;
+        bool final_check = false; 
+
+        for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
+            its_since_last_switch++;
+            if (its_since_last_switch >= min_its_after_switch){
+                const double wqTotal = this->primary_variables_.getWQTotal().value();
+                changed = this->updateWellControlAndStatusLocalIteration (ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger); 
+                if (changed){
+                    its_since_last_switch = 0;
+                    switch_count++;
+                }
+                if (!changed && final_check) {
+                    break;
+                } else {
+                    final_check = false;
+                }
+            }
+
+            assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+
+            const BVectorWell dx_well = this->linSys_.solve();
+
+            if (it > this->param_.strict_inner_iter_wells_) {
+                relax_convergence = true;
+                this->regularize_ = true;
+            }
+
+            const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
+            converged = report.converged();
+            if (converged) {
+                // if equations are sufficiently linear they might converge in less than min_its_after_switch
+                // in this case, make sure all constraints are satisfied before returning
+                if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
+                    final_check = true;
+                    its_since_last_switch = min_its_after_switch;
+                } else {
+                    break;
+                }
+            }
+
+            // getFinteWellResiduals returns false for nan/inf residuals
+            {
+                const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+                if (!isFinite)
+                    return false;
+
+                residual_history.push_back(residuals);
+            }
+
+            if (!converged) {
+                measure_history.push_back(this->getResidualMeasureValue(well_state,
+                                                                        residual_history[it],
+                                                                        this->param_.tolerance_wells_,
+                                                                        this->param_.tolerance_pressure_ms_wells_,
+                                                                        deferred_logger));
+
+                bool is_oscillate = false;
+                bool is_stagnate = false;
+
+                this->detectOscillations(measure_history, is_oscillate, is_stagnate);
+                // TODO: maybe we should have more sophisticated strategy to recover the relaxation factor,
+                // for example, to recover it to be bigger
+
+                if (is_oscillate || is_stagnate) {
+                    // HACK!
+                    std::string message;
+                    if (relaxation_factor == min_relaxation_factor) {
+                        ++stagnate_count;
+                        if (false) { // this disables the usage of the relaxed tolerance
+                            fmt::format_to(std::back_inserter(message), " Well {} observes severe stagnation and/or oscillation."
+                                                                        " We relax the tolerance and check for convergence. \n", this->name());
+                            const auto reportStag = getWellConvergence(summary_state, well_state, Base::B_avg_,
+                                                                       deferred_logger, true);
+                            if (reportStag.converged()) {
+                                converged = true;
+                                fmt::format_to(std::back_inserter(message), " Well {}  manages to get converged with relaxed tolerances in {} inner iterations", this->name(), it);
+                                deferred_logger.debug(message);
+                                return converged;
+                            }
+                        }
+                    }
+
+                    // a factor value to reduce the relaxation_factor
+                    constexpr double reduction_mutliplier = 0.9;
+                    relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
+
+                    // debug output
+                    if (is_stagnate) {
+                        fmt::format_to(std::back_inserter(message), " well {} observes stagnation in inner iteration {}\n", this->name(), it);
+                    }
+                    if (is_oscillate) {
+                        fmt::format_to(std::back_inserter(message), " well {} observes oscillation in inner iteration {}\n", this->name(), it);
+                    }
+                    fmt::format_to(std::back_inserter(message), " relaxation_factor is {} now\n", relaxation_factor);
+
+                    this->regularize_ = true;
+                    deferred_logger.debug(message);
+                }
+            }
+            updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
+            initPrimaryVariablesEvaluation();
+        }
+
+        if (converged) {
+            if (allow_switching){
+                // update operability if status change
+                const bool is_stopped = this->wellIsStopped();
+                if (this->wellHasTHPConstraints(summary_state)){
+                    this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
+                    this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
+                } else {
+                    this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
+                }
+                // We reset the well status to it's original state. Status is updated 
+                // on the outside based on operability status
+                this->wellStatus_ = well_status;
+            }
+            std::string message = fmt::format("   Well {} converged in {} inner iterations ("
+                                                    "{} control/status switches).", this->name(), it, switch_count);
+            if (relax_convergence) {
+                message.append(fmt::format("   (A relaxed tolerance was used after {} iterations)",
+                                           this->param_.strict_inner_iter_wells_));
+            }
+            deferred_logger.debug(message);
+        } else {
+            const std::string message = fmt::format("   Well {} did not converged in {} inner iterations ("
+                                                    "{} control/status switches).", this->name(), it, switch_count);
+            deferred_logger.debug(message);
+        }
+
+        return converged;
+    }
 
 
     template<typename TypeTag>

opm/simulators/wells/StandardWell.hpp

@@ -205,6 +205,15 @@ namespace Opm
                                       const GroupState& group_state,
                                       DeferredLogger& deferred_logger) override;
 
+        // iterate well equations including control switching
+        bool iterateWellEqWithSwitching(const Simulator& ebosSimulator,
+                                        const double dt,
+                                        const Well::InjectionControls& inj_controls,
+                                        const Well::ProductionControls& prod_controls,
+                                        WellState& well_state,
+                                        const GroupState& group_state,
+                                        DeferredLogger& deferred_logger) override;
+
         /// \brief Wether the Jacobian will also have well contributions in it.
         virtual bool jacobianContainsWellContributions() const override
         {

opm/simulators/wells/StandardWell_impl.hpp

@@ -2155,6 +2155,103 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    iterateWellEqWithSwitching(const Simulator& ebosSimulator,
+                               const double dt,
+                               const Well::InjectionControls& inj_controls,
+                               const Well::ProductionControls& prod_controls,
+                               WellState& well_state,
+                               const GroupState& group_state,
+                               DeferredLogger& deferred_logger)
+    {
+        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 
+        constexpr int min_its_after_switch = 2;
+        int its_since_last_switch = min_its_after_switch;
+        int switch_count= 0;
+        const auto well_status = this->wellStatus_;
+        const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
+        bool changed = 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); 
+                if (changed){
+                    its_since_last_switch = 0;
+                    switch_count++;
+                }
+                if (!changed && final_check) {
+                    break;
+                } else {
+                    final_check = false;
+                }
+            }
+  
+            assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+
+            if (it > this->param_.strict_inner_iter_wells_) {
+                relax_convergence = true;
+                this->regularize_ = true;
+            }
+
+            auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
+
+            converged = report.converged();
+            if (converged) {
+                // if equations are sufficiently linear they might converge in less than min_its_after_switch
+                // in this case, make sure all constraints are satisfied before returning
+                if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
+                    final_check = true;
+                    its_since_last_switch = min_its_after_switch;
+                } else {
+                    break;
+                }   
+            }
+
+            ++it;
+            solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
+            initPrimaryVariablesEvaluation();
+        } while (it < max_iter);
+        
+        if (converged) {
+            if (allow_switching){
+                // update operability if status change
+                const bool is_stopped = this->wellIsStopped();
+                if (this->wellHasTHPConstraints(summary_state)){
+                    this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
+                    this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
+                } else {
+                    this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
+                }
+                // We reset the well status to its original state. Status is updated
+                // on the outside based on operability status
+                // \Note for future reference: For the well to update its status to stop/shut,
+                // the flag changed_to_stopped_this_step_ in prepareWellBeforeAssembling needs to be set to true.
+                // 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;
+            }
+        } else {
+            const std::string message = fmt::format("   Well {} did not converged in {} inner iterations ("
+                                                    "{} control/status switches).", this->name(), it, switch_count);
+            deferred_logger.debug(message);
+            // add operability here as well ?
+        }
+        return converged;
+    }
+
     template<typename TypeTag>
     std::vector<double>
     StandardWell<TypeTag>::

opm/simulators/wells/WellConstraints.cpp

@@ -39,13 +39,16 @@ checkIndividualConstraints(SingleWellState& ws,
                            const SummaryState& summaryState,
                            const RateConvFunc& calcReservoirVoidageRates,
                            bool& thp_limit_violated_but_not_switched,
-                           DeferredLogger& deferred_logger) const
+                           DeferredLogger& deferred_logger,
+                           const std::optional<Well::InjectionControls>& inj_controls,
+                           const std::optional<Well::ProductionControls>& prod_controls) const
 {
     if (well_.isProducer()) {
         auto new_cmode = this->activeProductionConstraint(ws, summaryState,
                                                           calcReservoirVoidageRates,
                                                           thp_limit_violated_but_not_switched,
-                                                          deferred_logger);
+                                                          deferred_logger, 
+                                                          prod_controls);
         if (new_cmode != ws.production_cmode) {
             ws.production_cmode = new_cmode;
             return true;
@@ -55,7 +58,8 @@ checkIndividualConstraints(SingleWellState& ws,
     if (well_.isInjector()) {
         auto new_cmode = this->activeInjectionConstraint(ws, summaryState,
                                                         thp_limit_violated_but_not_switched,
-                                                        deferred_logger);
+                                                        deferred_logger, 
+                                                        inj_controls);
         if (new_cmode != ws.injection_cmode) {
             ws.injection_cmode = new_cmode;
             return true;
@@ -69,11 +73,12 @@ Well::InjectorCMode WellConstraints::
 activeInjectionConstraint(const SingleWellState& ws,
                           const SummaryState& summaryState,
                           bool& thp_limit_violated_but_not_switched,
-                          DeferredLogger& deferred_logger) const
+                          DeferredLogger& deferred_logger,
+                          const std::optional<Well::InjectionControls>& inj_controls) const
 {
     const PhaseUsage& pu = well_.phaseUsage();
 
-    const auto controls = well_.wellEcl().injectionControls(summaryState);
+    const auto controls =  inj_controls.has_value() ? inj_controls.value() : well_.wellEcl().injectionControls(summaryState);
     const auto currentControl = ws.injection_cmode;
 
     if (controls.hasControl(Well::InjectorCMode::BHP) && currentControl != Well::InjectorCMode::BHP)
@@ -166,10 +171,11 @@ activeProductionConstraint(const SingleWellState& ws,
                            const SummaryState& summaryState,
                            const RateConvFunc& calcReservoirVoidageRates,
                            bool& thp_limit_violated_but_not_switched,
-                           DeferredLogger& deferred_logger) const
+                           DeferredLogger& deferred_logger,
+                           const std::optional<Well::ProductionControls>& prod_controls) const
 {
     const PhaseUsage& pu = well_.phaseUsage();
-    const auto controls = well_.wellEcl().productionControls(summaryState);
+    const auto controls = prod_controls.has_value() ? prod_controls.value() : well_.wellEcl().productionControls(summaryState);
     const auto currentControl = ws.production_cmode;
 
     if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP) {

opm/simulators/wells/WellConstraints.hpp

@@ -24,9 +24,12 @@
 #ifndef OPM_WELL_CONSTRAINTS_HEADER_INCLUDED
 #define OPM_WELL_CONSTRAINTS_HEADER_INCLUDED
 
+#include <opm/input/eclipse/Schedule/Well/Well.hpp>
+
 #include <functional>
 #include <utility>
 #include <vector>
+#include <optional>
 
 namespace Opm
 {
@@ -55,21 +58,25 @@ public:
                                const SummaryState& summaryState,
                                const RateConvFunc& calcReservoirVoidageRates,
                                bool& thp_limit_violated_but_not_switched,
-                               DeferredLogger& deferred_logger) const;
+                               DeferredLogger& deferred_logger,
+                               const std::optional<Well::InjectionControls>& inj_controls = std::nullopt,
+                               const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
 
 private:
     WellInjectorCMode
     activeInjectionConstraint(const SingleWellState& ws,
                               const SummaryState& summaryState,
                               bool& thp_limit_violated_but_not_switched,
-                              DeferredLogger& deferred_logger) const;
+                              DeferredLogger& deferred_logger,
+                              const std::optional<Well::InjectionControls>& inj_controls = std::nullopt) const;
 
     WellProducerCMode
     activeProductionConstraint(const SingleWellState& ws,
                                const SummaryState& summaryState,
                                const RateConvFunc& calcReservoirVoidageRates,
                                bool& thp_limit_violated_but_not_switched,
-                               DeferredLogger& deferred_logger) const;
+                               DeferredLogger& deferred_logger,
+                               const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
 
     const WellInterfaceGeneric& well_; //!< Reference to well interface
 };

opm/simulators/wells/WellInterface.hpp

@@ -241,6 +241,14 @@ public:
                            const GroupState& group_state,
                            DeferredLogger& deferred_logger) /* const */;
 
+    bool updateWellControlAndStatusLocalIteration(const Simulator& ebos_simulator,
+                                                  WellState& well_state,
+                                                  const GroupState& group_state,
+                                                  const Well::InjectionControls& inj_controls,
+                                                  const Well::ProductionControls& prod_controls,
+                                                  const double WQTotal,
+                                                  DeferredLogger& deferred_logger);
+
     virtual void updatePrimaryVariables(const SummaryState& summary_state,
                                         const WellState& well_state,
                                         DeferredLogger& deferred_logger) = 0;
@@ -302,6 +310,9 @@ public:
                                const WellState& well_state,
                                DeferredLogger& deferred_logger);
 
+    bool updateWellOperabilityFromWellEq(const Simulator& ebos_simulator,
+                                         const WellState& well_state,
+                                         DeferredLogger& deferred_logger);
 
     // update perforation water throughput based on solved water rate
     virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
@@ -390,6 +401,14 @@ protected:
                                           const GroupState& group_state,
                                           DeferredLogger& deferred_logger) = 0;
 
+    virtual bool iterateWellEqWithSwitching(const Simulator& ebosSimulator,
+                                            const double dt,
+                                            const WellInjectionControls& inj_controls,
+                                            const WellProductionControls& prod_controls,
+                                            WellState& well_state,
+                                            const GroupState& group_state,
+                                            DeferredLogger& deferred_logger) = 0;
+
     bool iterateWellEquations(const Simulator& ebosSimulator,
                               const double dt,
                               WellState& well_state,

opm/simulators/wells/WellInterfaceFluidSystem.cpp

@@ -105,7 +105,9 @@ bool
 WellInterfaceFluidSystem<FluidSystem>::
 checkIndividualConstraints(SingleWellState& ws,
                            const SummaryState& summaryState,
-                           DeferredLogger& deferred_logger) const
+                           DeferredLogger& deferred_logger,
+                           const std::optional<Well::InjectionControls>& inj_controls,
+                           const std::optional<Well::ProductionControls>& prod_controls) const
 {
     auto rRates = [this](const int fipreg,
                          const int pvtRegion,
@@ -119,7 +121,7 @@ checkIndividualConstraints(SingleWellState& ws,
     return WellConstraints(*this).
             checkIndividualConstraints(ws, summaryState, rRates,
                                        this->operability_status_.thp_limit_violated_but_not_switched,
-                                       deferred_logger);
+                                       deferred_logger, inj_controls, prod_controls);
 }
 
 template <typename FluidSystem>

opm/simulators/wells/WellInterfaceFluidSystem.hpp

@@ -79,7 +79,9 @@ protected:
 
     bool checkIndividualConstraints(SingleWellState& ws,
                                     const SummaryState& summaryState,
-                                    DeferredLogger& deferred_logger) const;
+                                    DeferredLogger& deferred_logger,
+                                    const std::optional<Well::InjectionControls>& inj_controls = std::nullopt,
+                                    const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
 
     bool checkGroupConstraints(WellState& well_state,
                                const GroupState& group_state,

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -344,12 +344,25 @@ void WellInterfaceGeneric::setVFPProperties(const VFPProperties* vfp_properties_
     vfp_properties_ = vfp_properties_arg;
 }
 
-void WellInterfaceGeneric::setPrevSurfaceRates(WellState& well_state, 
+void WellInterfaceGeneric::setPrevSurfaceRates(WellState& well_state,
                                                const WellState& prev_well_state) const
     {
         auto& ws = well_state.well(this->index_of_well_);
-        if (!this->changedToOpenThisStep()){
+        auto& ws_prev = prev_well_state.well(this->index_of_well_);
-            ws.prev_surface_rates = prev_well_state.well(this->index_of_well_).surface_rates;
+        // The logic here is a bit fragile:
+        // We need non-zero prev_surface_rates for the purpose of providing explicit fractions
+        // (if needed) for vfp interpolation.
+        // We assume that current surface rates either are initialized from previous step
+        // or (if newly opened) from updateWellStateRates. This is fine unless well was
+        // stopped in previous step in which case it's rates will be zero. In this case,
+        // we select the previous rates of the previous well state (and hope for the best).
+        const bool zero_rates = std::all_of(ws.surface_rates.begin(), ws.surface_rates.end(),
+                [](double rate) {
+                    return rate == 0.; // TODO: should we use a threshhold for comparison?
+                } );
+
+        if (zero_rates) {
+            ws.prev_surface_rates = ws_prev.prev_surface_rates;
         } else {
             ws.prev_surface_rates = ws.surface_rates;
         }
@@ -552,6 +565,11 @@ bool WellInterfaceGeneric::stopppedOrZeroRateTarget(const SummaryState& summary_
 
 }
 
+void WellInterfaceGeneric::resetWellOperability()
+{
+    this->operability_status_.resetOperability();
+}
+
 double WellInterfaceGeneric::wmicrobes_() const
 {
     auto injectorType = this->well_ecl_.injectorType();

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -218,6 +218,8 @@ public:
         inj_fc_multiplier_ = inj_fc_multiplier;
     }
 
+    void resetWellOperability();
+
 protected:
     bool getAllowCrossFlow() const;
 

opm/simulators/wells/WellInterface_impl.hpp

@@ -254,7 +254,82 @@ namespace Opm
         return changed;
     }
 
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    updateWellControlAndStatusLocalIteration(const Simulator& ebos_simulator,
+                                             WellState& well_state,
+                                             const GroupState& group_state,
+                                             const Well::InjectionControls& inj_controls,
+                                             const Well::ProductionControls& prod_controls,
+                                             const double wqTotal,
+                                             DeferredLogger& deferred_logger)
+    {
+        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;
+        }
 
+        const double sgn = this->isInjector() ? 1.0 : -1.0;
+        if (!this->wellIsStopped()){
+            if (wqTotal*sgn <= 0.0){
+                this->stopWell();
+                return true;
+            } else {
+                bool changed = false;
+                auto& ws = well_state.well(this->index_of_well_);
+                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);
+                // TODO: with current way, the checkGroupConstraints might overwrite the result from checkIndividualConstraints, which remains to be investigated
+                if (hasGroupControl) {
+                    changed = this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
+                }
+
+                if (changed) {
+                    const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
+                                                                     ws.production_cmode == Well::ProducerCMode::THP;
+                    if (!thp_controlled){
+                        // don't call for thp since this might trigger additional local solve
+                        updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
+                    } else {
+                        ws.thp = this->getTHPConstraint(summary_state);
+                    }  
+                    updatePrimaryVariables(summary_state, well_state, deferred_logger);
+                }
+                return changed;
+            }
+        } else {
+            // well is stopped, check if current bhp allows reopening
+            const double bhp = well_state.well(this->index_of_well_).bhp;
+            double prod_limit = prod_controls.bhp_limit;
+            double inj_limit = inj_controls.bhp_limit;
+            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 
+                // option is probably minimum bhp on current vfp-curve, but some more functionality 
+                // 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);
+                if (this->isInjector()){
+                    inj_limit = std::min(bhp_thp, inj_controls.bhp_limit);
+                } else {
+                    prod_limit = std::max(bhp_thp, prod_controls.bhp_limit);
+                }
+            }
+            const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
+            if (bhp_diff > 0){
+                this->openWell();
+                return true;
+            } else {
+                return false;
+            }
+        }
+    }
 
     template<typename TypeTag>
     void
@@ -363,7 +438,13 @@ namespace Opm
         const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
         bool converged = false;
         try {
-            converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+            // TODO: the following two functions will be refactored to be one to reduce the code duplication
+            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);
+            }
+                 
         } 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);
@@ -545,13 +626,21 @@ namespace Opm
         const bool well_operable = this->operability_status_.isOperableAndSolvable();
 
         if (!well_operable && old_well_operable) {
-            if (this->well_ecl_.getAutomaticShutIn()) {
+            if (this->param_.local_well_solver_control_switching_) {
-                deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
+                deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
+                this->stopWell();
+                changed_to_stopped_this_step_ = true;
             } else {
-                if (!this->wellIsStopped()) {
+                // \Note: keep the old manner for now for testing checking.
-                    deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
+                // Will be investgiated and fixed in a later PR
-                    this->stopWell();
+                if (this->well_ecl_.getAutomaticShutIn()) {
-                    changed_to_stopped_this_step_ = true;
+                    deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
+                } else {
+                    if (!this->wellIsStopped()) {
+                        deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
+                        this->stopWell();
+                        changed_to_stopped_this_step_ = true;
+                    }
                 }
             }
         } else if (well_operable && !old_well_operable) {
@@ -706,7 +795,16 @@ 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) {
+                return;
+            } else {
+                deferred_logger.debug("Operability check using well equations did not converge for well "
+                                      + this->name() + ", reverting to classical approach." );
+            }
+        }
         this->operability_status_.resetOperability();
 
         bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
@@ -727,6 +825,23 @@ namespace Opm
         }
     }
 
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    updateWellOperabilityFromWellEq(const Simulator& ebos_simulator,
+                                    const WellState& well_state,
+                                    DeferredLogger& deferred_logger)
+    {
+        // only makes sense if we're using this parameter is true
+        assert(this->param_.local_well_solver_control_switching_);
+        this->operability_status_.resetOperability();
+        WellState well_state_copy = well_state;
+        const auto& group_state = ebos_simulator.problem().wellModel().groupState();
+        const double dt = ebos_simulator.timeStepSize();
+        // 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