Maximum allowed oscillation for groups part 2

Set a parameter for maximum allowed oscillation of group controls
The parameter is only used for iter > nupcol

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -72,6 +72,7 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     check_well_operability_ = Parameters::Get<Parameters::EnableWellOperabilityCheck>();
     check_well_operability_iter_ = Parameters::Get<Parameters::EnableWellOperabilityCheckIter>();
     max_number_of_well_switches_ = Parameters::Get<Parameters::MaximumNumberOfWellSwitches>();
+    max_number_of_group_switches_ = Parameters::Get<Parameters::MaximumNumberOfGroupSwitches>();
     use_average_density_ms_wells_ = Parameters::Get<Parameters::UseAverageDensityMsWells>();
     local_well_solver_control_switching_ = Parameters::Get<Parameters::LocalWellSolveControlSwitching>();
     use_implicit_ipr_ = Parameters::Get<Parameters::UseImplicitIpr>();
@@ -205,6 +206,8 @@ void BlackoilModelParameters<Scalar>::registerParameters()
         ("Enable the well operability checking during iterations");
     Parameters::Register<Parameters::MaximumNumberOfWellSwitches>
         ("Maximum number of times a well can switch to the same control");
+    Parameters::Register<Parameters::MaximumNumberOfGroupSwitches>
+        ("Maximum number of times a group can switch to the same control");
     Parameters::Register<Parameters::UseAverageDensityMsWells>
         ("Approximate segment densitities by averaging over segment and its outlet");
     Parameters::Register<Parameters::LocalWellSolveControlSwitching>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -121,6 +121,7 @@ template<class Scalar>
 struct RelaxedPressureTolMsw { static constexpr Scalar value = 1e4; };
 
 struct MaximumNumberOfWellSwitches { static constexpr int value = 3; };
+struct MaximumNumberOfGroupSwitches { static constexpr int value = 3; };
 struct UseAverageDensityMsWells { static constexpr bool value = false; };
 struct LocalWellSolveControlSwitching { static constexpr bool value = true; };
 struct UseImplicitIpr { static constexpr bool value = true; };
@@ -272,9 +273,12 @@ public:
     /// Whether to check well operability during iterations
     bool check_well_operability_iter_;
 
-    /// Maximum number of times a well can switch to the same controt
+    /// Maximum number of times a well can switch to the same control
     int max_number_of_well_switches_;
 
+    /// Maximum number of times group can switch to the same control
+    int max_number_of_group_switches_;
+
     /// Whether to approximate segment densities by averaging over segment and its outlet
     bool use_average_density_ms_wells_;
 

opm/simulators/wells/BlackoilWellModelConstraints.cpp

@@ -561,7 +561,8 @@ template<class Scalar>
 bool BlackoilWellModelConstraints<Scalar>::
 updateGroupIndividualControl(const Group& group,
                              const int reportStepIdx,
-                             std::map<std::pair<std::string,Phase>,std::string>& switched_inj,
+                             const int max_number_of_group_switch,
+                             std::map<std::pair<std::string,Phase>,std::vector<std::string>>& switched_inj,
                              std::map<std::string, std::vector<std::string>>& switched_prod,
                              std::map<std::string, std::pair<std::string, std::string>>& closed_offending_wells,
                              GroupState<Scalar>& group_state,
@@ -576,13 +577,43 @@ updateGroupIndividualControl(const Group& group,
             if (!group.hasInjectionControl(phase)) {
                 continue;
             }
+
+            bool group_is_oscillating = false;
+            if (switched_inj.count({group.name(), phase}) > 0) {
+                for (const auto& key : switched_inj[{group.name(), phase}]) {
+                    if (std::count(switched_inj[{group.name(), phase}].begin(), switched_inj[{group.name(), phase}].end(), key) >= max_number_of_group_switch) {
+                        const size_t sz = switched_inj[{group.name(), phase}].size();
+                        // only output if the two last keys are not the same.
+                        if (key !=  switched_inj[{group.name(), phase}][sz-2]) {
+                            if (wellModel_.comm().rank() == 0) {
+                                std::ostringstream os;
+                                os << phase;
+                                const std::string msg =
+                                    fmt::format("Group control for {} injector group {} is oscillating. Group control kept at {}.",
+                                                std::move(os).str(),
+                                                group.name(),
+                                                key);
+                                deferred_logger.info(msg);
+                            }
+                            switched_inj[{group.name(), phase}].push_back(key);
+                        }
+                        group_is_oscillating = true;
+                        break;
+                    }
+                }
+            }
+            if (group_is_oscillating) {
+                continue;
+            }
+
             const auto& changed_this = this->checkGroupInjectionConstraints(group,
                                                                             reportStepIdx,
                                                                             phase);
             if (changed_this.first != Group::InjectionCMode::NONE)
             {
-                switched_inj.insert_or_assign({group.name(), phase},
+                switched_inj[{group.name(), phase}].push_back(
-                                     Group::InjectionCMode2String(changed_this.first));
+                             Group::InjectionCMode2String(changed_this.first));
+
                 this->actionOnBrokenConstraints(group, changed_this.first, phase,
                                                 group_state, deferred_logger);
                 WellGroupHelpers<Scalar>::updateWellRatesFromGroupTargetScale(changed_this.second,
@@ -597,6 +628,27 @@ updateGroupIndividualControl(const Group& group,
         }
     }
     if (group.isProductionGroup()) {
+
+        if (switched_prod.count(group.name()) > 0) {
+            for (const auto& key : switched_prod[group.name()]) {
+                if (std::count(switched_prod[group.name()].begin(), switched_prod[group.name()].end(), key) >= max_number_of_group_switch) {
+                    const size_t sz = switched_prod[group.name()].size();
+                    // only output on rank 0 and if the two last keys are not the same.
+                    if (key !=  switched_prod[group.name()][sz-2]) {
+                        if (wellModel_.comm().rank() == 0) {
+                            const std::string msg =
+                            fmt::format("Group control for production group {} is oscillating. Group control kept at {}.",
+                                        group.name(),
+                                        key);
+                            deferred_logger.info(msg);
+                        }
+                        switched_prod[group.name()].push_back(key);
+                    }
+                    return false;
+                }
+            }
+        }
+
         const auto& changed_this = this->checkGroupProductionConstraints(group,
                                                                          reportStepIdx,
                                                                          deferred_logger);

opm/simulators/wells/BlackoilWellModelConstraints.hpp

@@ -73,7 +73,8 @@ public:
     //! \brief Update the individual controls for wells in a group. Return true if a group control is changed
     bool updateGroupIndividualControl(const Group& group,
                                       const int reportStepIdx,
-                                      std::map<std::pair<std::string,Phase>,std::string>& switched_inj,
+                                      const int max_number_of_group_switch,
+                                      std::map<std::pair<std::string,Phase>,std::vector<std::string>>& switched_inj,
                                       std::map<std::string, std::vector<std::string>>& switched_prod,
                                       std::map<std::string, std::pair<std::string, std::string>>& closed_offending_wells,
                                       GroupState<Scalar>& group_state,

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -574,7 +574,8 @@ template<class Scalar>
 bool BlackoilWellModelGeneric<Scalar>::
 checkGroupHigherConstraints(const Group& group,
                             DeferredLogger& deferred_logger,
-                            const int reportStepIdx)
+                            const int reportStepIdx,
+                            const int max_number_of_group_switch)
 {
     // Set up coefficients for RESV <-> surface rate conversion.
     // Use the pvtRegionIdx from the top cell of the first well.
@@ -626,6 +627,34 @@ checkGroupHigherConstraints(const Group& group,
         }
         const Phase all[] = { Phase::WATER, Phase::OIL, Phase::GAS };
         for (Phase phase : all) {
+            bool group_is_oscillating = false;
+            if (switched_inj_groups_.count({group.name(), phase}) > 0) {
+                for (const auto& key : switched_inj_groups_[{group.name(), phase}]) {
+                    if (std::count(switched_inj_groups_[{group.name(), phase}].begin(), switched_inj_groups_[{group.name(), phase}].end(), key) >= max_number_of_group_switch) {
+                        const size_t sz = switched_inj_groups_[{group.name(), phase}].size();
+                        // only output if the two last keys are not the same.
+                        if (key !=  switched_inj_groups_[{group.name(), phase}][sz-2]) {
+                            if (comm_.rank() == 0 ) {
+                                std::ostringstream os;
+                                os << phase;
+                                const std::string msg =
+                                    fmt::format("Group control for {} injector group {} is oscillating. Group control kept at {}.",
+                                                std::move(os).str(),
+                                                group.name(),
+                                                key);
+                                deferred_logger.info(msg);
+                            }
+                            switched_inj_groups_[{group.name(), phase}].push_back(key);
+                        }
+                        group_is_oscillating = true;
+                        break;
+                    }
+                }
+            }
+            if (group_is_oscillating) {
+                continue;
+            }
+
             // Check higher up only if under individual (not FLD) control.
             auto currentControl = this->groupState().injection_control(group.name(), phase);
             if (currentControl != Group::InjectionCMode::FLD && group.injectionGroupControlAvailable(phase)) {
@@ -647,7 +676,7 @@ checkGroupHigherConstraints(const Group& group,
                                                                        resv_coeff_inj,
                                                                        deferred_logger);
                 if (is_changed) {
-                    switched_inj_groups_.insert_or_assign({group.name(), phase}, Group::InjectionCMode2String(Group::InjectionCMode::FLD));
+                    switched_inj_groups_[{group.name(), phase}].push_back(Group::InjectionCMode2String(Group::InjectionCMode::FLD));
                     BlackoilWellModelConstraints(*this).
                         actionOnBrokenConstraints(group, Group::InjectionCMode::FLD,
                                                   phase, this->groupState(),
@@ -671,6 +700,25 @@ checkGroupHigherConstraints(const Group& group,
         // So when checking constraints, current groups rate must also be subtracted it's reduction rate
         const std::vector<Scalar> reduction_rates = this->groupState().production_reduction_rates(group.name());
 
+        if (switched_prod_groups_.count(group.name()) > 0) {
+            for (const auto& key : switched_prod_groups_[group.name()]) {
+                if (std::count(switched_prod_groups_[group.name()].begin(), switched_prod_groups_[group.name()].end(), key) >= max_number_of_group_switch) {
+                    const size_t sz = switched_prod_groups_[group.name()].size();
+                    // only output on rank 0 and if the two last keys are not the same.
+                    if (key !=  switched_prod_groups_[group.name()][sz-2]) {
+                        if (comm_.rank() == 0) {
+                            const std::string msg =
+                            fmt::format("Group control for production group {} is oscillating. Group control kept at {}.",
+                                        group.name(),
+                                        key);
+                            deferred_logger.info(msg);
+                        }
+                        switched_prod_groups_[group.name()].push_back(key);
+                    }
+                    return false;
+                }
+            }
+        }
         for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
             const Scalar local_current_rate = WellGroupHelpers<Scalar>::sumWellSurfaceRates(group,
                                                                                             schedule(),

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -375,7 +375,8 @@ protected:
 
     bool checkGroupHigherConstraints(const Group& group,
                                      DeferredLogger& deferred_logger,
-                                     const int reportStepIdx);
+                                     const int reportStepIdx,
+                                     const int max_number_of_group_switch);
 
     void updateAndCommunicateGroupData(const int reportStepIdx,
                                        const int iterationIdx,
@@ -599,7 +600,7 @@ protected:
 
     // Store maps of group name and new group controls for output
     std::map<std::string, std::vector<std::string>> switched_prod_groups_;
-    std::map<std::pair<std::string, Phase>, std::string> switched_inj_groups_;
+    std::map<std::pair<std::string, Phase>, std::vector<std::string>> switched_inj_groups_;
     // Store map of group name and close offending well for output
     std::map<std::string, std::pair<std::string, std::string>> closed_offending_wells_;
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -766,11 +766,11 @@ namespace Opm {
 
                 const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
                 std::string from = Group::InjectionCMode2String(oldControl);
-                if (to != from) {
+                if (to.back() != from) {
                     std::string msg = "    Injection Group " + name
                     + " control mode changed from ";
                     msg += from;
-                    msg += " to " + to;
+                    msg += " to " + to.back();
                     local_deferredLogger.info(msg);
                 }
             }
@@ -2466,25 +2466,11 @@ namespace Opm {
                         const int reportStepIdx,
                         const int iterationIdx)
     {
-
-        if (this->switched_prod_groups_.count(group.name()) > 0) {
-            for (const auto& key : this->switched_prod_groups_[group.name()]) {
-                if (std::count(this->switched_prod_groups_[group.name()].begin(), this->switched_prod_groups_[group.name()].end(), key) >= 3) {
-                    if (std::count(this->switched_prod_groups_[group.name()].begin(), this->switched_prod_groups_[group.name()].end(), key) == 3) {
-                        const std::string msg =
-                            fmt::format("Group control for group {} is oscilating. Group control kept at {}.",
-                                        group.name(),
-                                        key);
-                        deferred_logger.info(msg);
-                        this->switched_prod_groups_[group.name()].push_back(key);
-                    }
-                    return false;
-                }
-            }
-        }
-
         bool changed = false;
-        bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx);
+        // restrict the number of group switches but only after nupcol iterations.
+        const int nupcol = this->schedule()[reportStepIdx].nupcol();
+        const int max_number_of_group_switches = iterationIdx <= nupcol ? 9999 : param_.max_number_of_group_switches_;
+        bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx, max_number_of_group_switches);
         if (changed_hc) {
             changed = true;
             updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
@@ -2494,6 +2480,7 @@ namespace Opm {
             BlackoilWellModelConstraints(*this).
                 updateGroupIndividualControl(group,
                                              reportStepIdx,
+                                             max_number_of_group_switches,
                                              this->switched_inj_groups_,
                                              this->switched_prod_groups_,
                                              this->closed_offending_wells_,

tests/test_RestartSerialization.cpp

@@ -333,7 +333,7 @@ public:
         nupcol_wgstate_ = WGState<double>::serializationTestObject(dummy);
         last_glift_opt_time_ = 5.0;
         switched_prod_groups_ = {{"test4", {"test5", "test6"}}};
-        switched_inj_groups_ = {{{"test4", Phase::SOLVENT}, "test5"}};
+        switched_inj_groups_ = {{{"test4", Phase::SOLVENT}, {"test5", "test6"}}};
         closed_offending_wells_ = {{"test4", {"test5", "test6"}}};
     }