Merge pull request #5753 from totto82/maxGroupChange

Alternative way of avoiding group control oscillations

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,8 +561,9 @@ 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::string, std::string>& switched_prod,
+                             std::map<std::string, std::array<std::vector<Group::InjectionCMode>, 3>>& switched_inj,
+                             std::map<std::string, std::vector<Group::ProductionCMode>>& switched_prod,
                              std::map<std::string, std::pair<std::string, std::string>>& closed_offending_wells,
                              GroupState<Scalar>& group_state,
                              WellState<Scalar>& well_state,
@@ -576,13 +577,44 @@ updateGroupIndividualControl(const Group& group,
             if (!group.hasInjectionControl(phase)) {
                 continue;
             }
+            bool group_is_oscillating = false;
+            if (auto groupPos = switched_inj.find(group.name()); groupPos != switched_inj.end()) {
+                auto& ctrls = groupPos->second[static_cast<std::underlying_type_t<Phase>>(phase)];
+                for (const auto& ctrl : ctrls) {
+                    if (std::count(ctrls.begin(), ctrls.end(), ctrl) < max_number_of_group_switch) {
+                        continue;
+                    }
+
+                    if (ctrls.back() != *(ctrls.end() - 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(),
+                                            Group::InjectionCMode2String(ctrl));
+                            deferred_logger.info(msg);
+                        }
+                        ctrls.push_back(ctrl);
+                    }
+                    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()][static_cast<std::underlying_type_t<Phase>>(phase)].push_back(
-                                     Group::InjectionCMode2String(changed_this.first));
+                             changed_this.first);
+
                 this->actionOnBrokenConstraints(group, changed_this.first, phase,
                                                 group_state, deferred_logger);
                 WellGroupHelpers<Scalar>::updateWellRatesFromGroupTargetScale(changed_this.second,
@@ -597,6 +629,28 @@ updateGroupIndividualControl(const Group& group,
         }
     }
     if (group.isProductionGroup()) {
+
+        if (auto groupPos = switched_prod.find(group.name()); groupPos != switched_prod.end()) {
+            auto& ctrls = groupPos->second;
+            for (const auto& ctrl : ctrls) {
+                if (std::count(ctrls.begin(), ctrls.end(), ctrl) < max_number_of_group_switch) {
+                    continue;
+                }
+
+                if (ctrls.back() != *(ctrls.end() - 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(),
+                                    Group::ProductionCMode2String(ctrl));
+                        deferred_logger.info(msg);
+                    }
+                    ctrls.push_back(ctrl);
+                }
+                return false;
+            }
+        }
+
         const auto& changed_this = this->checkGroupProductionConstraints(group,
                                                                          reportStepIdx,
                                                                          deferred_logger);
@@ -612,8 +666,7 @@ updateGroupIndividualControl(const Group& group,
                                             group_state, deferred_logger);
 
             if(changed) {
-                switched_prod.insert_or_assign(group.name(),
+                switched_prod[group.name()].push_back(changed_this.first);
-                                    Group::ProductionCMode2String(changed_this.first));
                 WellGroupHelpers<Scalar>::updateWellRatesFromGroupTargetScale(changed_this.second,
                                                                               group,
                                                                               wellModel_.schedule(),

opm/simulators/wells/BlackoilWellModelConstraints.hpp

@@ -73,8 +73,9 @@ 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::string, std::string>& switched_prod,
+                                      std::map<std::string, std::array<std::vector<Group::InjectionCMode>, 3>>& switched_inj,
+                                      std::map<std::string, std::vector<Group::ProductionCMode>>& switched_prod,
                                       std::map<std::string, std::pair<std::string, std::string>>& closed_offending_wells,
                                       GroupState<Scalar>& group_state,
                                       WellState<Scalar>& well_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,36 @@ checkGroupHigherConstraints(const Group& group,
         }
         const Phase all[] = { Phase::WATER, Phase::OIL, Phase::GAS };
         for (Phase phase : all) {
+            bool group_is_oscillating = false;
+            if (auto groupPos = switched_inj_groups_.find(group.name()); groupPos != switched_inj_groups_.end()) {
+                auto& ctrls = groupPos->second[static_cast<std::underlying_type_t<Phase>>(phase)];
+                for (const auto& ctrl : ctrls) {
+                    if (std::count(ctrls.begin(), ctrls.end(), ctrl) < max_number_of_group_switch) {
+                        continue;
+                    }
+
+                    if (ctrls.back() != *(ctrls.end() - 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(),
+                                            Group::InjectionCMode2String(ctrl));
+                            deferred_logger.info(msg);
+                        }
+                        ctrls.push_back(ctrl);
+                    }
+                    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 +678,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()][static_cast<std::underlying_type_t<Phase>>(phase)].push_back(Group::InjectionCMode::FLD);
                     BlackoilWellModelConstraints(*this).
                         actionOnBrokenConstraints(group, Group::InjectionCMode::FLD,
                                                   phase, this->groupState(),
@@ -671,6 +702,26 @@ 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 (auto groupPos = switched_prod_groups_.find(group.name()); groupPos != switched_prod_groups_.end()) {
+            auto& ctrls = groupPos->second;
+            for (const auto& ctrl : ctrls) {
+                if (std::count(ctrls.begin(), ctrls.end(), ctrl) < max_number_of_group_switch) {
+                    continue;
+                }
+
+                if (ctrls.back() != *(ctrls.end() - 2)) {
+                    if (comm_.rank() == 0) {
+                        const std::string msg =
+                        fmt::format("Group control for production group {} is oscillating. Group control kept at {}.",
+                                    group.name(),
+                                    Group::ProductionCMode2String(ctrl));
+                        deferred_logger.info(msg);
+                    }
+                    ctrls.push_back(ctrl);
+                }
+                return false;
+            }
+        }
         for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
             const Scalar local_current_rate = WellGroupHelpers<Scalar>::sumWellSurfaceRates(group,
                                                                                             schedule(),
@@ -714,8 +765,7 @@ checkGroupHigherConstraints(const Group& group,
                                                   deferred_logger);
 
                 if (changed) {
-                    switched_prod_groups_.insert_or_assign(group.name(),
+                    switched_prod_groups_[group.name()].push_back(Group::ProductionCMode::FLD);
-                                                           Group::ProductionCMode2String(Group::ProductionCMode::FLD));
                     WellGroupHelpers<Scalar>::updateWellRatesFromGroupTargetScale(scaling_factor,
                                                                                   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,
@@ -598,8 +599,8 @@ protected:
     bool wellStructureChangedDynamically_{false};
 
     // Store maps of group name and new group controls for output
-    std::map<std::string, std::string> switched_prod_groups_;
+    std::map<std::string, std::vector<Group::ProductionCMode>> switched_prod_groups_;
-    std::map<std::pair<std::string, Phase>, std::string> switched_inj_groups_;
+    std::map<std::string, std::array<std::vector<Group::InjectionCMode>, 3>> 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

@@ -749,10 +749,11 @@ namespace Opm {
         // report group switching
         if (this->terminal_output_) {
 
-            for (const auto& [name, to] : this->switched_prod_groups_) {
+            for (const auto& [name, ctrls] : this->switched_prod_groups_) {
                 const Group::ProductionCMode& oldControl = this->prevWGState().group_state.production_control(name);
                 std::string from = Group::ProductionCMode2String(oldControl);
-                if (to != from) {
+                std::string to = Group::ProductionCMode2String(ctrls.back());
+                if (ctrls.back() != oldControl) {
                     std::string msg = "    Production Group " + name
                     + " control mode changed from ";
                     msg += from;
@@ -760,18 +761,30 @@ namespace Opm {
                     local_deferredLogger.info(msg);
                 }
             }
-            for (const auto& [key, to] : this->switched_inj_groups_) {
-                const std::string& name = key.first;
-                const Opm::Phase& phase = key.second;
 
-                const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
+            for (const auto& [grname, grdata] : this->switched_inj_groups_) {
-                std::string from = Group::InjectionCMode2String(oldControl);
+                //const std::string& name = key.first;
-                if (to != from) {
+                //const Opm::Phase& phase = key.second;
-                    std::string msg = "    Injection Group " + name
+                const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS};
-                    + " control mode changed from ";
+                for (Phase phase : all) {
-                    msg += from;
+                    const auto& ctrls = grdata[static_cast<std::underlying_type_t<Phase>>(phase)];
-                    msg += " to " + to;
+
-                    local_deferredLogger.info(msg);
+                    if (ctrls.empty()) {
+                        continue;
+                    }
+                    if ( !this->prevWGState().group_state.has_injection_control(grname, phase))
+                        continue;
+
+                    const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(grname, phase);
+                    std::string from = Group::InjectionCMode2String(oldControl);
+                    std::string to = Group::InjectionCMode2String(ctrls.back());
+                    if (ctrls.back() != oldControl) {
+                        std::string msg = "    Injection Group " + grname
+                        + " control mode changed from ";
+                        msg += from;
+                        msg += " to " + to;
+                        local_deferredLogger.info(msg);
+                    }
                 }
             }
         }
@@ -2208,13 +2221,9 @@ namespace Opm {
 
         bool changed_well_group = false;
         // Check group individual constraints.
-        const int nupcol = this->schedule()[episodeIdx].nupcol();
+        const Group& fieldGroup = this->schedule().getGroup("FIELD", episodeIdx);
-        // don't switch group control when iterationIdx > nupcol
+        changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
-        // to avoid oscilations between group controls
+
-        if (iterationIdx <= nupcol) {
-            const Group& fieldGroup = this->schedule().getGroup("FIELD", episodeIdx);
-            changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
-        }
         // Check wells' group constraints and communicate.
         bool changed_well_to_group = false;
         {
@@ -2262,7 +2271,6 @@ namespace Opm {
         }
 
         // update wsolvent fraction for REIN wells
-        const Group& fieldGroup = this->schedule().getGroup("FIELD", episodeIdx);
         this->updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
 
         return { changed_well_group, more_network_update };
@@ -2472,7 +2480,10 @@ namespace Opm {
                         const int iterationIdx)
     {
         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);
@@ -2482,6 +2493,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_,

opm/simulators/wells/WellInterface_impl.hpp

@@ -213,8 +213,10 @@ namespace Opm
             from = WellProducerCMode2String(ws.production_cmode);
         }
         bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
-
+        const int episodeIdx = simulator.episodeIndex();
-        if (oscillating) {
+        const int iterationIdx = simulator.model().newtonMethod().numIterations();
+        const int nupcol = schedule[episodeIdx].nupcol();
+        if (oscillating && iterationIdx > nupcol) {
             // only output frist time
             bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
             if (output) {

tests/test_RestartSerialization.cpp

@@ -325,8 +325,9 @@ public:
         last_valid_wgstate_ = WGState<double>::serializationTestObject(dummy);
         nupcol_wgstate_ = WGState<double>::serializationTestObject(dummy);
         last_glift_opt_time_ = 5.0;
-        switched_prod_groups_ = {{"test4", "test5"}};
+        switched_prod_groups_ = {{"test4", {Group::ProductionCMode::NONE, Group::ProductionCMode::ORAT}}};
-        switched_inj_groups_ = {{{"test4", Phase::SOLVENT}, "test5"}};
+        const auto controls = {Group::InjectionCMode::NONE, Group::InjectionCMode::RATE, Group::InjectionCMode::RATE };
+        switched_inj_groups_ = {{"test4", {controls, {}, controls} }};
         closed_offending_wells_ = {{"test4", {"test5", "test6"}}};
     }