From afdfe55bbc98196bb3261f7f82bc3f3d532097a4 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Thu, 15 Apr 2021 08:14:52 +0200
Subject: [PATCH] restrict maximum number of times a well can switch to the
 same control

---
 .../flow/BlackoilModelParametersEbos.hpp      | 21 +++++++++++++++----
 opm/simulators/wells/WellInterface.hpp        |  2 ++
 opm/simulators/wells/WellInterface_impl.hpp   | 19 +++++++++++++++++
 3 files changed, 38 insertions(+), 4 deletions(-)

diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
index d1eafdc04..4a8631e2d 100644
--- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
@@ -160,6 +160,10 @@ template<class TypeTag, class MyTypeTag>
 struct AlternativeWellRateInit {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct MaximumNumberOfWellSwitches {
+    using type = UndefinedProperty;
+};
 
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
@@ -300,6 +304,10 @@ struct RelaxedPressureTolMsw<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 0.5e5;
 };
+template<class TypeTag>
+struct MaximumNumberOfWellSwitches<TypeTag, TTag::FlowModelParameters> {
+    static constexpr int value = 3;
+};
 
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
@@ -402,14 +410,18 @@ namespace Opm
         /// The file name of the deck
         std::string deck_file_name_;
 
-        // Whether to add influences of wells between cells to the matrix and preconditioner matrix
+        /// Whether to add influences of wells between cells to the matrix and preconditioner matrix
         bool matrix_add_well_contributions_;
 
-        // Whether to check well operability
+        /// Whether to check well operability
         bool check_well_operability_;
-        // Whether to check well operability during iterations
+        /// Whether to check well operability during iterations
         bool check_well_operability_iter_;
 
+        /// Maximum number of times a well can switch to the same controt
+        int max_number_of_well_switches_;
+
+
 
         /// Construct from user parameters or defaults.
         BlackoilModelParametersEbos()
@@ -444,7 +456,7 @@ namespace Opm
             matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
             check_well_operability_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
             check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
-
+            max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
             deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
         }
 
@@ -482,6 +494,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices");
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking");
             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");
         }
     };
 } // namespace Opm
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 28a7042b6..ddc4229d6 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -278,6 +278,8 @@ protected:
 
     bool changed_to_stopped_this_step_ = false;
 
+    std::vector< std::string> well_control_log_;
+
     double wpolymer() const;
 
     double wfoam() const;
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 8159d4220..4c72f9c55 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -61,6 +61,7 @@ namespace Opm
                 }
             }
         }
+        well_control_log_.clear();
     }
 
 
@@ -172,7 +173,24 @@ namespace Opm
         } else {
             from = Well::ProducerCMode2String(ws.production_cmode);
         }
+        bool oscillating = std::count(well_control_log_.begin(), well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
 
+        if (oscillating) {
+            // only output frist time
+            bool output = std::count(well_control_log_.begin(), well_control_log_.end(), from) == param_.max_number_of_well_switches_;
+            if (output) {
+                std::ostringstream ss;
+                ss << "    The control model for well " << this->name()
+                   << " is oscillating\n"
+                   << "    We don't allow for more than "
+                   << param_.max_number_of_well_switches_
+                   << " switches. The control is kept at " << from;
+                deferred_logger.info(ss.str());
+                // add one more to avoid outputting the same info again
+                well_control_log_.push_back(from);
+            }
+            return false;
+        }
         bool changed = false;
         if (iog == IndividualOrGroup::Individual) {
             changed = this->checkIndividualConstraints(ws, summaryState);
@@ -193,6 +211,7 @@ namespace Opm
             } else {
                 to = Well::ProducerCMode2String(ws.production_cmode);
             }
+            well_control_log_.push_back(from);
             std::ostringstream ss;
             ss << "    Switching control mode for well " << this->name()
                << " from " << from