Add parameter for restricting status oscillations between stop and open

The parameters can improve convergence for cases where wells jums between stop and open

Default is 99 i.e. by default we dont restrict ossiclations

We have two levels of parameters for this
The first is due restict the number in the inner iterations
(This functionality is already in the code)
--max-well-status-switch-in-inner-iter-wells

The second is to restrict the number during a time-step
--max-well-status-switch-iter-wells

The second could lead to wells being shut too early and should be
used more carefully

flowexperimental/flowexp.hpp

@@ -169,6 +169,7 @@ public:
         Parameters::Hide<Parameters::MaxNewtonIterationsWithInnerWellIterations>();
         Parameters::Hide<Parameters::MaxInnerIterWells>();
         Parameters::Hide<Parameters::MaxWellStatusSwitchInInnerIterWells>();
+        Parameters::Hide<Parameters::MaxWellStatusSwitchForWells>();
         Parameters::Hide<Parameters::MaxSinglePrecisionDays<Scalar>>();
         Parameters::Hide<Parameters::MinStrictCnvIter>();
         Parameters::Hide<Parameters::MinStrictMbIter>();

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -62,7 +62,8 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     max_niter_inner_well_iter_ = Parameters::Get<Parameters::MaxNewtonIterationsWithInnerWellIterations>();
     shut_unsolvable_wells_ = Parameters::Get<Parameters::ShutUnsolvableWells>();
     max_inner_iter_wells_ = Parameters::Get<Parameters::MaxInnerIterWells>();
-    max_well_status_switch_ = Parameters::Get<Parameters::MaxWellStatusSwitchInInnerIterWells>();
+    max_well_status_switch_inner_iter_ = Parameters::Get<Parameters::MaxWellStatusSwitchInInnerIterWells>();
+    max_well_status_switch_ = Parameters::Get<Parameters::MaxWellStatusSwitchForWells>();
     maxSinglePrecisionTimeStep_ = Parameters::Get<Parameters::MaxSinglePrecisionDays<Scalar>>() * 24 * 60 * 60;
     min_strict_cnv_iter_ = Parameters::Get<Parameters::MinStrictCnvIter>();
     min_strict_mb_iter_ = Parameters::Get<Parameters::MinStrictMbIter>();
@@ -185,7 +186,9 @@ void BlackoilModelParameters<Scalar>::registerParameters()
     Parameters::Register<Parameters::MaxInnerIterWells>
         ("Maximum number of inner iterations for standard wells");
     Parameters::Register<Parameters::MaxWellStatusSwitchInInnerIterWells>
-        ("Maximum number of status switching (shut<->open) in inner iterations for wells");
+        ("Maximum number of status switching (stop<->open) for a well during inner iterations.");
+    Parameters::Register<Parameters::MaxWellStatusSwitchForWells>
+        ("Maximum number of status switching (stop<->open) for a well during a time-step.");
     Parameters::Register<Parameters::AlternativeWellRateInit>
         ("Use alternative well rate initialization procedure");
     Parameters::Register<Parameters::RegularizationFactorWells<Scalar>>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -104,6 +104,7 @@ struct MaxNewtonIterationsWithInnerWellIterations { static constexpr int value =
 struct MaxInnerIterMsWells { static constexpr int value = 100; };
 struct MaxInnerIterWells { static constexpr int value = 50; };
 struct MaxWellStatusSwitchInInnerIterWells { static constexpr int value = 99; };
+struct MaxWellStatusSwitchForWells { static constexpr int value = 99; };
 struct ShutUnsolvableWells { static constexpr bool value = true; };
 struct AlternativeWellRateInit { static constexpr bool value = true; };
 struct StrictOuterIterWells { static constexpr int value = 6; };
@@ -324,7 +325,10 @@ public:
     /// Maximum number of iterations in the well/group switch algorithm
     int well_group_constraints_max_iterations_;
 
-    /// Maximum number of status switches (open<->shut> in local well iterations
+    /// Maximum number of status switches (open<->stop> in local well iterations
+    int max_well_status_switch_inner_iter_;
+
+    /// Maximum number of status switches (open<->stop> during a time step
     int max_well_status_switch_;
 
     /// Nonlinear solver type: newton or nldd.

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1634,7 +1634,7 @@ namespace Opm
         // that unless this number is >1, we may get stuck in a cycle 
         const int min_its_after_switch = 3;
         // We also want to restrict the number of status switches to avoid oscillation between STOP<->OPEN
-        const int max_status_switch = this->param_.max_well_status_switch_;
+        const int max_status_switch = this->param_.max_well_status_switch_inner_iter_;
         int its_since_last_switch = min_its_after_switch;
         int switch_count= 0;
         int status_switch_count = 0;

opm/simulators/wells/StandardWell_impl.hpp

@@ -2386,7 +2386,7 @@ namespace Opm
         // that unless this number is >1, we may get stuck in a cycle
         constexpr int min_its_after_switch = 4;
         // We also want to restrict the number of status switches to avoid oscillation between STOP<->OPEN
-        const int max_status_switch = this->param_.max_well_status_switch_;
+        const int max_status_switch = this->param_.max_well_status_switch_inner_iter_;
         int its_since_last_switch = min_its_after_switch;
         int switch_count= 0;
         // if we fail to solve eqs, we reset status/operability before leaving

opm/simulators/wells/WellInterface.hpp

@@ -180,7 +180,6 @@ public:
                                      const GroupState<Scalar>& group_state,
                                      DeferredLogger& deferred_logger);
 
-
     virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
                                          const Scalar& bhp,
                                          std::vector<Scalar>& well_flux,
@@ -368,6 +367,7 @@ protected:
     std::vector<Scalar> B_avg_;
     bool changed_to_stopped_this_step_ = false;
     bool thp_update_iterations = false;
+    int number_of_well_reopenings_{0};
 
     Scalar wpolymer() const;
     Scalar wfoam() const;

opm/simulators/wells/WellInterface_impl.hpp

@@ -574,7 +574,6 @@ namespace Opm
     {
         OPM_TIMEFUNCTION();
         const auto& summary_state = simulator.vanguard().summaryState();
-        bool is_operable = true;
         bool converged = true;
         auto& ws = well_state.well(this->index_of_well_);
         // if well is stopped, check if we can reopen
@@ -585,10 +584,13 @@ namespace Opm
                 // no intersection with ipr
                 const auto msg = fmt::format("estimateOperableBhp: Did not find operable BHP for well {}", this->name());
                 deferred_logger.debug(msg);
-                is_operable = false;
+                // well can't operate using explicit fractions stop the well
                 // solve with zero rates
-                solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
+                converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
                 this->stopWell();
+                this->operability_status_.can_obtain_bhp_with_thp_limit = false;
+                this->operability_status_.obey_thp_limit_under_bhp_limit = false;
+                return converged;
             } else {
                 // solve well with the estimated target bhp (or limit)
                 ws.thp = this->getTHPConstraint(summary_state);
@@ -598,9 +600,8 @@ namespace Opm
             }
         }
         // solve well-equation
-        if (is_operable) {
-            converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
-        }
+        converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+
 
         const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
         // check stability of solution under thp-control
@@ -633,11 +634,13 @@ namespace Opm
             this->openWell();
             auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
             if (!bhp_target.has_value()) {
-                // well can't operate using explicit fractions
-                is_operable = false;
                 // solve with zero rate
+                // well can't operate using explicit fractions stop the well
                 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
                 this->stopWell();
+                this->operability_status_.can_obtain_bhp_with_thp_limit = false;
+                this->operability_status_.obey_thp_limit_under_bhp_limit = false;
+                return converged;
             } else {
                 // solve well with the estimated target bhp (or limit)
                 const Scalar bhp = std::max(bhp_target.value(),
@@ -655,9 +658,8 @@ namespace Opm
             }
         }
         // update operability
-        is_operable = is_operable && !this->wellIsStopped();
-        this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
-        this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
+        this->operability_status_.can_obtain_bhp_with_thp_limit = !this->wellIsStopped();
+        this->operability_status_.obey_thp_limit_under_bhp_limit = !this->wellIsStopped();
         return converged;
     }
 
@@ -897,6 +899,25 @@ namespace Opm
                             [](Scalar rate) { return rate != Scalar(0.0); });
 
             this->operability_status_.solvable = true;
+            if (number_of_well_reopenings_ >= this->param_.max_well_status_switch_) {
+                // only output the first time
+                if (number_of_well_reopenings_ == this->param_.max_well_status_switch_) {
+                    const std::string msg = fmt::format("well {} is oscillating between open and stop. \n"
+                                                    "We don't allow for more than {} re-openings "
+                                                    "and the well is therefore kept stopped.",
+                                                     this->name(), number_of_well_reopenings_);
+                    deferred_logger.debug(msg);
+                }
+                this->stopWell();
+                changed_to_stopped_this_step_ = true;
+                bool converged_zero_rate = this->solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
+                if (this->param_.shut_unsolvable_wells_ && !converged_zero_rate ) {
+                    this->operability_status_.solvable = false;
+                }
+                // we increse the number of reopenings to avoid output in the next iteration
+                number_of_well_reopenings_++;
+                return;
+            }
             bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
 
             if (converged) {
@@ -908,6 +929,7 @@ namespace Opm
                     this->operability_status_.resetOperability();
                     this->openWell();
                     deferred_logger.debug("    " + this->name() + " is re-opened after being stopped during local solve");
+                    number_of_well_reopenings_++;
                 }
             } else {
                 // unsolvable wells are treated as not operable and will not be solved for in this iteration.