updates according to Atgeirrs comments

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1312,8 +1312,8 @@ namespace Opm
     {
         // 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)
+        // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target)
-        // where Eq(x)=0 is the well equation setup with bhp control and primary varables x 
+        // where Eq(x)=0 is the well equation setup with bhp control and primary variables x 
 
         // We shouldn't have zero rates at this stage, but check
         bool zero_rates;
@@ -1363,6 +1363,7 @@ namespace Opm
             const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
             const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
             for (size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
+                // well primary variable derivatives in EvalWell start at position Indices::numEq 
                 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
             }
             ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
@@ -1610,7 +1611,9 @@ namespace Opm
         this->regularize_ = false;
         const auto& summary_state = ebosSimulator.vanguard().summaryState();
 
-        // Max status switch frequency should be 2 to avoid getting stuck in cycle
+        // Always take a few (more than one) iterations after a switch before allowing a new switch
+        // The optimal number here is subject to further investigation, but it has been observerved 
+        // that unless this number is >1, we may get stuck in a cycle 
         const int min_its_after_switch = 3;
         int its_since_last_switch = min_its_after_switch;
         int switch_count= 0;
@@ -1754,7 +1757,7 @@ namespace Opm
         } else {
             this->wellStatus_ = well_status_orig;
             this->operability_status_ = operability_orig;            
-            const std::string message = fmt::format("   Well {} did not converged in {} inner iterations ("
+            const std::string message = fmt::format("   Well {} did not converge in {} inner iterations ("
                                                     "{} control/status switches).", this->name(), it, switch_count);
             deferred_logger.debug(message);
         }

opm/simulators/wells/StandardWell.hpp

@@ -242,7 +242,7 @@ namespace Opm
             const double alq_value,
             DeferredLogger& deferred_logger) const override;
 
-        void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger);                
+        void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger) override;
 
         virtual void computeWellRatesWithBhp(
             const Simulator& ebosSimulator,

opm/simulators/wells/StandardWell_impl.hpp

@@ -841,8 +841,8 @@ namespace Opm
     {   
         // 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)
+        // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target)
-        // where Eq(x)=0 is the well equation setup with bhp control and primary varables x 
+        // where Eq(x)=0 is the well equation setup with bhp control and primary variables x 
 
         // We shouldn't have zero rates at this stage, but check
         bool zero_rates;
@@ -853,7 +853,7 @@ namespace Opm
         }
         auto& ws = well_state.well(this->index_of_well_);
         if (zero_rates) {
-            const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be probelmatic", this->name());
+            const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
             deferred_logger.debug(msg);
             /*
             // could revert to standard approach here:    
@@ -882,7 +882,7 @@ namespace Opm
         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();
+        const size_t nEq = this->primary_variables_.numWellEq();
         BVectorWell rhs(1);
         rhs[0].resize(nEq);
         // rhs = 0 except -1 for control eq
@@ -899,6 +899,7 @@ namespace Opm
             EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
             const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
             for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
+                // well primary variable derivatives in EvalWell start at position Indices::numEq 
                 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
             }
             ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
@@ -2321,14 +2322,15 @@ namespace Opm
                                const bool fixed_status /*false*/)
     {
         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 
+        // Always take a few (more than one) iterations after a switch before allowing a new switch
+        // The optimal number here is subject to further investigation, but it has been observerved 
+        // that unless this number is >1, we may get stuck in a cycle 
         constexpr int min_its_after_switch = 4;
         int its_since_last_switch = min_its_after_switch;
         int switch_count= 0;
@@ -2341,7 +2343,7 @@ namespace Opm
         bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
         allow_switching = allow_switching && (!fixed_control || !fixed_status);
         bool changed = false;
-        bool final_check = false; 
+        bool final_check = false;
         // well needs to be set operable or else solving/updating of re-opened wells is skipped
         this->operability_status_.resetOperability();
         this->operability_status_.solvable = true;
@@ -2364,7 +2366,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_) {
@@ -2383,7 +2385,7 @@ namespace Opm
                     its_since_last_switch = min_its_after_switch;
                 } else {
                     break;
-                }   
+                }
             }
 
             ++it;
@@ -2411,7 +2413,7 @@ namespace Opm
         } else {
             this->wellStatus_ = well_status_orig;
             this->operability_status_ = operability_orig;
-            const std::string message = fmt::format("   Well {} did not converged in {} inner iterations ("
+            const std::string message = fmt::format("   Well {} did not converge in {} inner iterations ("
                                                     "{} 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

@@ -543,7 +543,7 @@ intersectWithIPR(const VFPProdTable& table,
 {   
     // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection 
     // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable 
-    // intersections, the one corresonding the largest flo i returned.
+    // intersections, the one corresponding the largest flo is returned.
     // The adjust_bhp-function is used to adjust the vfp-table bhp-values to actual bhp-values due
     // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. 
 

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -962,7 +962,7 @@ getFloIPR(const WellState& well_state,
     const double& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
     const double& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
     // The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type), 
-    // but we can equally use it to pick/add the corresonding ipr_a, ipr_b  
+    // but we can equally use it to pick/add the corresponding ipr_a, ipr_b  
     return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), 
                           detail::getFlo(table, aqua_b, liquid_b, vapour_b));
 }

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -748,7 +748,7 @@ void WellInterfaceGeneric::
 prepareForPotentialCalculations(const SummaryState& summary_state,
                                 WellState& well_state, 
                                 Well::InjectionControls& inj_controls,
-                                Well::ProductionControls& prod_controls)
+                                Well::ProductionControls& prod_controls) const
 {
     const bool has_thp = this->wellHasTHPConstraints(summary_state);
     auto& ws = well_state.well(this->index_of_well_);

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -248,7 +248,7 @@ protected:
     void prepareForPotentialCalculations(const SummaryState& summary_state,
                                          WellState& well_state, 
                                          Well::InjectionControls& inj_controls,
-                                         Well::ProductionControls& prod_controls);                                  
+                                         Well::ProductionControls& prod_controls) const;
 
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {

opm/simulators/wells/WellInterface_impl.hpp

@@ -270,7 +270,7 @@ namespace Opm
     {
         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;
         }
@@ -278,7 +278,6 @@ namespace Opm
         const double sgn = this->isInjector() ? 1.0 : -1.0;
         if (!this->wellIsStopped()){
             if (wqTotal*sgn <= 0.0 && !fixed_status){
-                //std::cout << "Stopping well:" << this->name() << std::endl;
                 this->stopWell();
                 return true;
             } else {
@@ -322,13 +321,10 @@ namespace Opm
                     // if the well can operate, it must at least be able to produce at the lowest bhp of the bhp-curve (explicit fractions)
                     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 = (this->isInjector())? inj_limit : prod_limit;
                 if (has_thp) {
@@ -357,9 +353,7 @@ 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();
@@ -367,7 +361,12 @@ namespace Opm
         initPrimaryVariablesEvaluation();
 
         if (this->isProducer()) {
-            gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
+            const auto& schedule = simulator.vanguard().schedule();
+            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;
@@ -462,7 +461,7 @@ namespace Opm
                     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);
@@ -512,32 +511,29 @@ namespace Opm
 
         const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
         // check stability of solution under thp-control
-        if (true) { 
+        if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && isThp) {
-            if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && isThp) {
+            auto rates = well_state.well(this->index_of_well_).surface_rates;
-                auto rates = well_state.well(this->index_of_well_).surface_rates;
+            this->adaptRatesForVFP(rates);
-                this->adaptRatesForVFP(rates);
+            this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger);
-                this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger);
+            bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state);
-                bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state);
+            if (!is_stable) {
-                if (!is_stable) {
+                // solution converged to an unstable point! 
-                    // solution converged to an unstable point! 
+                this->operability_status_.use_vfpexplicit = true;
-                    this->operability_status_.use_vfpexplicit = true;
+                auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
-                    // msg = ...
+                // if we find an intersection with a sufficiently lower bhp, re-solve equations
-                    auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
+                const double reltol = 1e-3;
-                    // if we find an intersection with a sufficiently lower bhp, re-solve equations
+                const double cur_bhp = ws.bhp;
-                    const double reltol = 1e-3;
+                if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
-                    const double cur_bhp = ws.bhp;
+                    const auto msg = fmt::format("Well {} converged to an unstable solution, re-solving", this->name());
-                    if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
+                    deferred_logger.debug(msg);
-                        const auto msg = fmt::format("isStableSolution: Well {} converged to an unstable solution, re-solving", this->name());
+                    solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger);
-                        deferred_logger.debug(msg);
+                    // re-solve with hopefully good initial guess
-                        solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger);
+                    ws.thp = this->getTHPConstraint(summary_state);
-                        // re-solve with hopefully good initial guess
+                    converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
-                        ws.thp = this->getTHPConstraint(summary_state);
+                }
-                        converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+            }
-                    }
-                }          
-            }  
         }
-        
+
         if (!converged) {
             // Well did not converge, switch to explicit fractions
             this->operability_status_.use_vfpexplicit = true;
@@ -900,8 +896,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.has_well(well_name)) {
+        if(glo.active() && glo.has_well(well_name)) {
-            auto increment = glo.gaslift_increment();
+            const auto increment = glo.gaslift_increment();
             auto alq = well_state.getALQ(well_name);
             bool converged;
             while (alq > 0) {
@@ -934,9 +930,8 @@ namespace Opm
             deferred_logger.info(msg);
             return;
         }
-        const auto& well_ecl = this->wellEcl();
         const auto& schedule = ebos_simulator.vanguard().schedule();
-        auto report_step_idx = ebos_simulator.episodeIndex();
+        const 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(
@@ -952,6 +947,7 @@ 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();
@@ -970,7 +966,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) {
@@ -1016,7 +1012,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