reworked implementation of group control feasibility

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -122,6 +122,7 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
 
     nupcol_group_rate_tolerance_ = Parameters::Get<Parameters::NupcolGroupRateTolerance<Scalar>>();
     well_group_constraints_max_iterations_ = Parameters::Get<Parameters::WellGroupConstraintsMaxIterations>();
+    group_control_fraction_tolerance_ = Parameters::Get<Parameters::GroupControlFractionTolerance<Scalar>>();
 }
 
 template<class Scalar>
@@ -314,8 +315,10 @@ void BlackoilModelParameters<Scalar>::registerParameters()
     Parameters::Hide<Parameters::DebugEmitCellPartition>();
 
     Parameters::Register<Parameters::WellGroupConstraintsMaxIterations>
-    ("Maximum number of iterations in the well/group switching algorithm");
+        ("Maximum number of iterations in the well/group switching algorithm");
 
+    Parameters::Register<Parameters::GroupControlFractionTolerance<Scalar>>
+        ("Tolerance for minimal allowed fraction for controlled phase before employing a fallback control");
     // if openMP is available, use two threads per mpi rank by default
 #if _OPENMP
     Parameters::SetDefault<Parameters::ThreadsPerProcess>(2);

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -159,6 +159,8 @@ struct NewtonMinIterations { static constexpr int value = 2; };
 
 struct WellGroupConstraintsMaxIterations { static constexpr int value = 1; };
 template<class Scalar>
+struct GroupControlFractionTolerance { static constexpr Scalar value = 1e-4; };
+template<class Scalar>
 struct LocalToleranceScalingMb { static constexpr Scalar value = 1.0; };
 
 template<class Scalar>
@@ -354,6 +356,9 @@ public:
     /// Maximum number of status switches (open<->stop> during a time step
     int max_well_status_switch_;
 
+    /// Lower tolerance for the guiderate ratio or well-fraction of active phase before switching to group fallback control
+    Scalar group_control_fraction_tolerance_;
+
     /// Nonlinear solver type: newton or nldd
     std::string nonlinear_solver_;
 

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1318,10 +1318,14 @@ updateAndCommunicateGroupData(const int reportStepIdx,
     this->wellState().communicateGroupRates(comm_);
     this->groupState().communicate_rates(comm_);
 
+    if (iterCtx.isFirstGlobalIteration()) {
+        group_state_helper.updatePreviousGroupProductionRates(fieldGroup);
+    }
+
     if (update_wellgrouptarget) {
         for (const auto& well : well_container_generic_) {
-            const auto& ws = this->wellState().well(well->indexOfWell());
-            const auto& group = this->schedule().getGroup(well->wellEcl().groupName(), well->currentStep());
+            auto& ws = this->wellState().well(well->indexOfWell());
+            const Group& group = this->schedule().getGroup(well->wellEcl().groupName(), reportStepIdx);
             std::vector<Scalar> resv_coeff(this->numPhases(), 0.0);
             const int fipnum = 0;
             const int pvtreg = well->pvtRegionIdx();
@@ -1332,7 +1336,10 @@ updateAndCommunicateGroupData(const int reportStepIdx,
             }
             const Scalar efficiencyFactor = well->wellEcl().getEfficiencyFactor() *
                                     ws.efficiency_scaling_factor;
-            auto& group_target = this->wellState().well(well->indexOfWell()).group_target;
+            auto& group_target = ws.group_target;
+            // ensure group target fallback and its flag are reset
+            ws.group_target_fallback = std::nullopt;
+            ws.use_group_target_fallback = false;
             if (well->isProducer()) {
                 group_target = group_state_helper.getWellGroupTargetProducer(
                     well->name(),
@@ -1346,7 +1353,26 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                     const std::string msg = fmt::format("Well {} is under GRUP control but no valid group target "
                         "could be determined. Switching the well to under BHP control.", well->name());
                     group_state_helper.deferredLogger().debug(msg);
-                    this->wellState().well(well->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
+                    ws.production_cmode = Well::ProducerCMode::BHP;
+                }
+                if (group_target.has_value()) {
+                    // check whether group giving target is operating at its original cmode, and if not, compute the target 
+                    // corresponding to original cmode for potential later use in case current cmode is/becomes infeasible
+                    const Group& target_group = this->schedule().getGroup(group_target->group_name, reportStepIdx);
+                    const auto cmode_orig = target_group.productionControls(summaryState_).cmode;
+                    if (cmode_orig != group_target->production_cmode &&
+                        cmode_orig != Group::ProductionCMode::FLD &&
+                        cmode_orig != Group::ProductionCMode::NONE) {
+                        auto& group_target_fallback = ws.group_target_fallback;
+                        group_target_fallback = group_state_helper.getWellGroupTargetProducer(
+                            well->name(),
+                            well->wellEcl().groupName(),
+                            group,
+                            ws.surface_rates.data(),
+                            efficiencyFactor,
+                            resv_coeff,
+                            cmode_orig);
+                    }
                 }
             } else {
                 const auto& well_controls = well->wellEcl().injectionControls(summaryState_);

opm/simulators/wells/FractionCalculator.hpp

@@ -54,6 +54,9 @@ public:
                     const bool always_include_this);
     Scalar localFraction(const std::string& name,
                          const std::string& always_included_child);
+    Scalar guideRate(const std::string& name,
+                     const std::string& always_included_child,
+                     const bool always_use_potentials);
 
 private:
     std::string parent(const std::string& name);
@@ -63,9 +66,6 @@ private:
     std::pair<Scalar,int> guideRateSum(const Group& group,
                         const std::string& always_included_child,
                         const bool always_use_potentials);
-    Scalar guideRate(const std::string& name,
-                     const std::string& always_included_child,
-                     const bool always_use_potentials);
     int groupControlledWells(const std::string& group_name,
                              const std::string& always_included_child);
     GuideRate::RateVector getGroupRateVector(const std::string& group_name);

opm/simulators/wells/GroupState.cpp

@@ -67,6 +67,7 @@ bool GroupState<Scalar>::operator==(const GroupState& other) const
            this->m_network_leaf_node_production_rates == other.m_network_leaf_node_production_rates &&
            this->production_controls == other.production_controls &&
            this->prod_red_rates == other.prod_red_rates &&
+           this->m_prev_production_rates == other.m_prev_production_rates &&
            this->inj_red_rates == other.inj_red_rates &&
            this->inj_resv_rates == other.inj_resv_rates &&
            this->inj_rein_rates == other.inj_rein_rates &&
@@ -120,6 +121,27 @@ GroupState<Scalar>::production_rates(const std::string& gname) const
     return group_iter->second;
 }
 
+template<class Scalar>
+const std::vector<Scalar>&
+GroupState<Scalar>::prev_production_rates(const std::string& gname) const
+{
+    auto group_iter = this->m_prev_production_rates.find(gname);
+    if (group_iter == this->m_prev_production_rates.end())
+        throw std::logic_error("No such group");
+
+    return group_iter->second;
+}
+
+template<class Scalar>
+void GroupState<Scalar>::update_prev_production_rates(const std::string& gname,
+                                                      const std::vector<Scalar>& rates)
+{
+    if (rates.size() != this->num_phases)
+        throw std::logic_error("Wrong number of phases");
+
+    this->m_prev_production_rates[gname] = rates;
+}
+
 template<class Scalar>
 const std::vector<Scalar>&
 GroupState<Scalar>::network_leaf_node_production_rates(const std::string& gname) const

opm/simulators/wells/GroupState.hpp

@@ -64,6 +64,9 @@ public:
                                            const std::vector<Scalar>& rates);
     const std::vector<Scalar>& production_reduction_rates(const std::string& gname) const;
 
+    void update_prev_production_rates(const std::string& gname, const std::vector<Scalar>& rates);
+    const std::vector<Scalar>& prev_production_rates(const std::string& gname) const;
+
     bool has_injection_reduction_rates(const std::string& gname) const;
     void update_injection_reduction_rates(const std::string& gname,
                                           const std::vector<Scalar>& rates);
@@ -213,6 +216,7 @@ public:
         serializer(m_production_rates);
         serializer(m_network_leaf_node_production_rates);
         serializer(production_controls);
+        serializer(m_prev_production_rates);
         serializer(group_thp);
         serializer(prod_red_rates);
         serializer(inj_red_rates);
@@ -234,6 +238,7 @@ private:
     std::map<std::string, std::vector<Scalar>> m_production_rates;
     std::map<std::string, std::vector<Scalar>> m_network_leaf_node_production_rates;
     std::map<std::string, Group::ProductionCMode> production_controls;
+    std::map<std::string, std::vector<Scalar>> m_prev_production_rates;
     std::map<std::string, std::vector<Scalar>> prod_red_rates;
     std::map<std::string, std::vector<Scalar>> inj_red_rates;
     std::map<std::string, std::vector<Scalar>> inj_surface_rates;

opm/simulators/wells/GroupStateHelper.cpp

@@ -540,6 +540,14 @@ GroupStateHelper<Scalar, IndexTraits>::
 getProductionGuideTargetMode(const Group& group) const
 {
     const auto cmode = this->groupState().production_control(group.name());
+    return this->getProductionGuideTargetModeFromControlMode(cmode);
+}
+
+template<typename Scalar, typename IndexTraits>
+GuideRateModel::Target
+GroupStateHelper<Scalar, IndexTraits>::
+getProductionGuideTargetModeFromControlMode(const Group::ProductionCMode cmode) const
+{
     switch (cmode) {
     case Group::ProductionCMode::ORAT:
         return GuideRateModel::Target::OIL;
@@ -767,17 +775,20 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::str
                                                              do_addback);
 
     // Avoid negative target rates coming from too large local reductions.
-    return GroupTarget{group.name(), std::max(Scalar(0.0), target / efficiency_factor)};
+    // guiderate_fraction is not used for injectors, so set to 0.0
+    return GroupTarget{group.name(), std::max(Scalar(0.0), target / efficiency_factor),
+                       Group::ProductionCMode::NONE, current_group_control, 0.0}; 
 }
 
 template <typename Scalar, typename IndexTraits>
 std::optional<typename SingleWellState<Scalar, IndexTraits>::GroupTarget>
 GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::string& name,
-                                                                 const std::string& parent,
-                                                                 const Group& group,
-                                                                 const Scalar* rates,
-                                                                 const Scalar efficiency_factor,
-                                                                 const std::vector<Scalar>& resv_coeff) const
+                                                                     const std::string& parent,
+                                                                     const Group& group,
+                                                                     const Scalar* rates,
+                                                                     const Scalar efficiency_factor,
+                                                                     const std::vector<Scalar>& resv_coeff,
+                                                                     std::optional<Group::ProductionCMode> fallback_cmode) const
 {
     // This function computes a wells group target.
     // 'parent' will be the name of 'group'. But if we recurse, 'name' and
@@ -787,10 +798,14 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
     // part of. Later it is the accumulated factor including the group efficiency factor
     // of the child of group.
     OPM_TIMEFUNCTION();
-    const Group::ProductionCMode& current_group_control = this->groupState().production_control(group.name());
+    const Group::ProductionCMode& current_group_cmode = this->groupState().production_control(group.name());
+    // If fallback_cmode is provided and different from current_group_cmode, we will translate the target from
+    // original mode to fallback mode using group fractions at previous step. Otherwise, fallback_cmode is 
+    // the same as current_group_cmode and no translation will be done.
+    const Group::ProductionCMode target_group_cmode = fallback_cmode.has_value() ? fallback_cmode.value() : current_group_cmode;
 
-    if (current_group_control == Group::ProductionCMode::FLD
-        || current_group_control == Group::ProductionCMode::NONE) {
+    if (current_group_cmode == Group::ProductionCMode::FLD
+        || current_group_cmode  == Group::ProductionCMode::NONE) {
         // Return if we are not available for parent group.
         if (!group.productionGroupControlAvailable()) {
             return std::nullopt;
@@ -802,7 +817,8 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
                                                 parent_group,
                                                 rates,
                                                 efficiency_factor * group.getGroupEfficiencyFactor(),
-                                                resv_coeff);
+                                                resv_coeff,
+                                                fallback_cmode);
     }
 
     // This can be false for FLD-controlled groups, we must therefore
@@ -813,17 +829,17 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
 
     // If we are here, we are at the topmost group to be visited in the recursion.
     // This is the group containing the control we will check against.
-
+    // Calculators for target mode
     GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{*this,
                                                                     resv_coeff,
-                                                                    group};
+                                                                    target_group_cmode};
 
     GroupStateHelpers::FractionCalculator<Scalar, IndexTraits> fcalc {this->schedule_,
                                                                       *this,
                                                                       this->summary_state_,
                                                                       this->report_step_,
                                                                       &this->guide_rate_,
-                                                                      this->getProductionGuideTargetMode(group),
+                                                                      this->getProductionGuideTargetModeFromControlMode(target_group_cmode),
                                                                       /*is_producer=*/true,
                                                                       /*injection_phase=*/Phase::OIL};
 
@@ -841,7 +857,20 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
         return fcalc.localFraction(child, always_included);
     };
 
-    const Scalar orig_target = this->getProductionGroupTarget(group);
+    Scalar orig_target = this->getProductionGroupTarget(group); // current cmode
+    if (target_group_cmode != current_group_cmode) {
+        // translate orig_target to target_cmode using "previous" group rate fractions
+        const auto& prev_rates = this->groupState().prev_production_rates(group.name());
+        // Calculator for current cmode:
+        GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc_current{*this,
+                                                                               resv_coeff,
+                                                                               current_group_cmode};
+        const Scalar prev_rate_current_cmode = tcalc_current.calcModeRateFromRates(prev_rates);
+        const Scalar prev_rate_target_cmode = tcalc.calcModeRateFromRates(prev_rates);
+        // shouldn't happen that prev_rate_current is zero since that is the controlled rate, but ...
+        const Scalar ratio = (prev_rate_current_cmode != 0.0) ? prev_rate_target_cmode / prev_rate_current_cmode : 1.0;
+        orig_target *= ratio; // now in target cmode
+    }
     // Switch sign since 'rates' are negative for producers.
     const Scalar current_rate_available = -tcalc.calcModeRateFromRates(rates);
     const auto chain = this->groupChainTopBot(name, group.name());
@@ -861,9 +890,18 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
                                                              local_reduction_lambda,
                                                              local_fraction_lambda,
                                                              do_addback);
-
+    // Finally, we provide the well-to-group guide-rate ratio for subsequent diagnostics of 
+    // target feasibility. Note that we cannot use the target directly for this, since a ~zero 
+    // target may very well just be the result of non-converged group tree. 
+    Scalar guide_rate_fraction = 1.0;
+    const Scalar well_guide_rate = fcalc.guideRate(name, name, /*always_use_potentials=*/false);
+    const Scalar group_guide_rate = fcalc.guideRate(group.name(), name, /*always_use_potentials=*/false);
+    if (group_guide_rate > 0.0) {
+        guide_rate_fraction = (well_guide_rate / group_guide_rate);
+    }  
     // Avoid negative target rates coming from too large local reductions.
-    return GroupTarget{group.name(), std::max(Scalar(0.0), target / efficiency_factor)};
+    return GroupTarget{group.name(), std::max(Scalar(0.0), target / efficiency_factor),
+                       target_group_cmode, Group::InjectionCMode::NONE, guide_rate_fraction};
 }
 
 template <typename Scalar, typename IndexTraits>
@@ -1198,6 +1236,19 @@ GroupStateHelper<Scalar, IndexTraits>::updateGroupProductionRates(const Group& g
     this->groupState().update_production_rates(group.name(), rates);
 }
 
+template <typename Scalar, typename IndexTraits>
+void
+GroupStateHelper<Scalar, IndexTraits>::updatePreviousGroupProductionRates(const Group& group)
+{
+    OPM_TIMEFUNCTION();
+    for (const std::string& group_name : group.groups()) {
+        const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_);
+        this->updatePreviousGroupProductionRates(group_tmp);
+    }
+    const std::vector<Scalar>& prev_rates = this->groupState().production_rates(group.name());
+    this->groupState().update_prev_production_rates(group.name(), prev_rates);
+}
+
 template <typename Scalar, typename IndexTraits>
 void
 GroupStateHelper<Scalar, IndexTraits>::updateGroupTargetReduction(const Group& group,

opm/simulators/wells/GroupStateHelper.hpp

@@ -51,7 +51,6 @@
 
 namespace Opm
 {
-
 template <typename Scalar, typename IndexTraits>
 class GroupStateHelper
 {
@@ -261,6 +260,9 @@ public:
     /// @return The GuideRateModel::Target based on the group's production control mode
     GuideRateModel::Target getProductionGuideTargetMode(const Group& group) const;
 
+    GuideRateModel::Target
+    getProductionGuideTargetModeFromControlMode(const Group::ProductionCMode cmode) const;
+
     std::pair<Scalar, Group::ProductionCMode>
     getAutoChokeGroupProductionTargetRate(const Group& bottom_group,
                                           const Group& group,
@@ -282,7 +284,8 @@ public:
                                                           const Group& group,
                                                           const Scalar* rates,
                                                           const Scalar efficiency_factor,
-                                                          const std::vector<Scalar>& resv_coeff) const;
+                                                          const std::vector<Scalar>& resv_coeff,
+                                                          std::optional<Group::ProductionCMode> target_cmode = std::nullopt) const;
 
     GuideRate::RateVector getWellRateVector(const std::string& name) const;
 
@@ -446,6 +449,8 @@ public:
 
     void updateGroupProductionRates(const Group& group);
 
+    void updatePreviousGroupProductionRates(const Group& group);
+
     void updateGroupTargetReduction(const Group& group,
                                     const bool is_injector);
 

opm/simulators/wells/MultisegmentWell.hpp

@@ -164,6 +164,9 @@ namespace Opm {
 
         int setPrimaryVars(typename std::vector<Scalar>::const_iterator it) override;
 
+        void getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
+                                    DeferredLogger& deferred_logger) const override;
+
     protected:
         // regularize msw equation
         bool regularize_;

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -702,6 +702,43 @@ outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const
     }
 }
 
+template<class FluidSystem, class Indices>
+void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
+scaledWellFractions(std::vector<Scalar>& fractions,
+                    DeferredLogger& /* deferred_logger */) const
+{
+    fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        fractions[oil_pos] = 1.0;
+        if (has_wfrac_variable) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            fractions[water_pos] = value_[0][WFrac];
+            fractions[oil_pos] -= fractions[water_pos];
+        }
+
+        if (has_gfrac_variable) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            fractions[gas_pos] = value_[0][GFrac];
+            fractions[oil_pos] -= fractions[gas_pos];
+        }
+    }
+    else if (has_wfrac_variable) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        fractions[water_pos] = 1.0;
+
+        if (has_gfrac_variable) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            fractions[gas_pos] = value_[0][GFrac];
+            fractions[water_pos] -= fractions[gas_pos];
+        }
+    }
+    else if (has_gfrac_variable) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        fractions[gas_pos] = 1.0;
+    }
+}
+
 #include <opm/simulators/utils/InstantiationIndicesMacros.hpp>
 
 INSTANTIATE_TYPE_INDICES(MultisegmentWellPrimaryVariables, double)

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -146,6 +146,10 @@ public:
     //! output the segments with pressure close to lower pressure limit for debugging purpose
     void outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const;
 
+    //! \brief get scaled fractions of top segment.
+    void scaledWellFractions(std::vector<Scalar>& fractions,
+                             DeferredLogger& deferred_logger) const;
+
 private:
     //! \brief Initialize evaluations from values.
     void setEvaluationsFromValues();

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1996,6 +1996,10 @@ namespace Opm
                     : groupStateHelper.groupState();
                 GroupStateHelperType groupStateHelper_copy = groupStateHelper;
                 auto group_guard = groupStateHelper_copy.pushGroupState(group_state);
+                // For production wells under group control, ensure feasibility before assembling control equation
+                if (this->wellUnderGroupControl(ws) && this->isProducer() && !stopped_or_zero_target) {
+                    this->updateGroupTargetFallbackFlag(well_state, deferred_logger);
+                }
                 MultisegmentWellAssemble(*this).
                         assembleControlEq(groupStateHelper_copy,
                                         inj_controls,
@@ -2351,6 +2355,16 @@ namespace Opm
         return num_seg * num_eq;
     }
 
+
+    template <typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
+                           DeferredLogger& deferred_logger) const
+    {
+        this->primary_variables_.scaledWellFractions(scaled_fractions, deferred_logger);
+    }
+
     template <typename TypeTag>
     typename MultisegmentWell<TypeTag>::FSInfo
     MultisegmentWell<TypeTag>::

opm/simulators/wells/SingleWellState.cpp

@@ -54,6 +54,7 @@ SingleWellState(const std::string& name_,
     , prev_surface_rates(pu.numActivePhases())
     , perf_data(perf_input.size(), !is_producer, pu.numActivePhases())
     , trivial_group_target(false)
+    , use_group_target_fallback(false)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
@@ -411,6 +412,8 @@ bool SingleWellState<Scalar, IndexTraits>::operator==(const SingleWellState& rhs
            this->alq_state == rhs.alq_state &&
            this->primaryvar == rhs.primaryvar &&
            this->group_target == rhs.group_target &&
+           this->group_target_fallback == rhs.group_target_fallback &&
+           this->use_group_target_fallback == rhs.use_group_target_fallback &&
            this->was_shut_before_action_applied == rhs.was_shut_before_action_applied;
 }
 

opm/simulators/wells/SingleWellState.hpp

@@ -25,6 +25,7 @@
 
 #include <opm/input/eclipse/Schedule/Well/WellEnums.hpp>
 #include <opm/input/eclipse/Schedule/Events.hpp>
+#include <opm/input/eclipse/Schedule/Group/Group.hpp>
 
 #include <opm/material/fluidsystems/PhaseUsageInfo.hpp>
 
@@ -86,6 +87,8 @@ public:
         serializer(primaryvar);
         serializer(alq_state);
         serializer(group_target);
+        serializer(group_target_fallback);
+        serializer(use_group_target_fallback);
         serializer(was_shut_before_action_applied);
     }
 
@@ -120,16 +123,26 @@ public:
 
     struct GroupTarget {
         std::string group_name;
-        Scalar target_value;
+        Scalar target_value{0.0};
+        Group::ProductionCMode production_cmode {Group::ProductionCMode::NONE};
+        Group::InjectionCMode injection_cmode {Group::InjectionCMode::NONE};
+        Scalar guiderate_ratio{1.0}; // well to group guide rate ratio for diagnostics
 
         bool operator==(const GroupTarget& other) const {
-            return group_name == other.group_name && target_value == other.target_value;
+            return (group_name == other.group_name 
+                 && target_value == other.target_value 
+                 && production_cmode == other.production_cmode
+                 && injection_cmode == other.injection_cmode
+                 && guiderate_ratio == other.guiderate_ratio);
         }
 
         template<class Serializer>
         void serializeOp(Serializer& serializer) {
             serializer(group_name);
             serializer(target_value);
+            serializer(production_cmode);
+            serializer(injection_cmode);
+            serializer(guiderate_ratio);
         }
     };
 
@@ -143,6 +156,8 @@ public:
     PerfData<Scalar> perf_data;
     bool trivial_group_target;
     std::optional<GroupTarget> group_target;
+    std::optional<GroupTarget> group_target_fallback;
+    bool use_group_target_fallback;
     SegmentState<Scalar> segments;
     Events events;
     WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};

opm/simulators/wells/StandardWell.hpp

@@ -249,6 +249,9 @@ namespace Opm
 
         int setPrimaryVars(typename std::vector<Scalar>::const_iterator it) override;
 
+        void getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
+                                    DeferredLogger& deferred_logger) const override;
+
     protected:
         bool regularize_;
 

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -692,6 +692,46 @@ processFractions()
     }
 }
 
+template<class FluidSystem, class Indices>
+void StandardWellPrimaryVariables<FluidSystem,Indices>::
+scaledWellFractions(std::vector<Scalar>& fractions,
+                    DeferredLogger& deferred_logger) const
+{
+    if constexpr (Indices::enableSolvent) {
+        OPM_DEFLOG_THROW(std::runtime_error,
+                         "Function scaledWellFractions does not support solvent yet.",
+                         deferred_logger);
+    }
+    fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        fractions[oil_pos] = 1.0;
+        if (has_wfrac_variable) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            fractions[water_pos] = value_[WFrac];
+            fractions[oil_pos] -= fractions[water_pos];
+        }
+        if (has_gfrac_variable) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            fractions[gas_pos] = value_[GFrac];
+            fractions[oil_pos] -= fractions[gas_pos];
+        }
+    }
+    else if (has_wfrac_variable) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        fractions[water_pos] = 1.0;
+        if (has_gfrac_variable) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            fractions[gas_pos] = value_[GFrac];
+            fractions[water_pos] -= fractions[gas_pos];
+        }
+    }
+    else if (has_gfrac_variable) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        fractions[gas_pos] = 1.0;
+    }
+}
+
 template<class FluidSystem, class Indices>
 typename StandardWellPrimaryVariables<FluidSystem,Indices>::Scalar
 StandardWellPrimaryVariables<FluidSystem,Indices>::

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -149,6 +149,10 @@ public:
     void setValue(const int idx, const Scalar val)
     { value_[idx] = val; }
 
+    //! \brief Get scaled well fraction values for all phases, based on the current primary variable values.
+    void scaledWellFractions(std::vector<Scalar>& fractions,
+                             DeferredLogger& deferred_logger) const;
+
 private:
     //! \brief Initialize evaluations from values.
     void setEvaluationsFromValues();

opm/simulators/wells/StandardWell_impl.hpp

@@ -476,6 +476,10 @@ namespace Opm
             auto& group_state = solving_with_zero_rate
                 ? empty_group_state
                 : groupStateHelper.groupState();
+            // For production wells under group control, ensure feasibility before assembling control equation
+            if (this->wellUnderGroupControl(ws) && this->isProducer() && !stopped_or_zero_target) {
+                this->updateGroupTargetFallbackFlag(well_state, deferred_logger);
+            }
             GroupStateHelperType groupStateHelper_copy = groupStateHelper;
             auto group_guard = groupStateHelper_copy.pushGroupState(group_state);
             StandardWellAssemble<FluidSystem,Indices>(*this).
@@ -2648,6 +2652,16 @@ namespace Opm
     }
 
 
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
+                           DeferredLogger& deferred_logger) const
+    {
+        this->primary_variables_.scaledWellFractions(scaled_fractions, deferred_logger);
+    }
+
+
     template <typename TypeTag>
     typename StandardWell<TypeTag>::Eval
     StandardWell<TypeTag>::

opm/simulators/wells/WellGroupControls.cpp

@@ -255,14 +255,15 @@ getGroupProductionControl(const Group& group,
     rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     // gconsale may adjust the grat target.
-    // the adjusted rates is send to the targetCalculator
 
-    GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{groupStateHelper, resv_coeff, group};
-
-    const auto target_rate = well_state.well(well_.indexOfWell()).group_target;
-    if (target_rate) {
+    const auto group_target = well_state.well(well_.indexOfWell()).use_group_target_fallback
+        ? well_state.well(well_.indexOfWell()).group_target_fallback
+        : well_state.well(well_.indexOfWell()).group_target;
+    if (group_target) {
+        // the adjusted rates is sent to the targetCalculator
+        GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{groupStateHelper, resv_coeff, group_target->production_cmode};
         const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
-        control_eq = current_rate - target_rate->target_value;
+        control_eq = current_rate - group_target->target_value;
     } else {
         const auto& controls = well.productionControls(summaryState);
         control_eq = bhp - controls.bhp_limit;
@@ -306,23 +307,25 @@ getGroupProductionTargetRate(const Group& group,
     std::vector<Scalar> resv_coeff(well_.phaseUsage().numActivePhases(), 1.0);
     rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
-    GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{groupStateHelper,
-                                                                   resv_coeff,
-                                                                   group};
+    const auto& ws = well_state.well(well_.indexOfWell());
+    const auto& group_target = ws.use_group_target_fallback
+        ? ws.group_target_fallback
+        : ws.group_target;
 
-    const auto target_rate = well_state.well(well_.indexOfWell()).group_target;
-    if (!target_rate) {
+    if (!group_target) {
         return 1.0;
     }
-    if (target_rate->target_value == 0.0) {
+    if (group_target->target_value == 0.0) {
         return 0.0;
     }
-    const auto& ws = well_state.well(well_.indexOfWell());
+    GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{groupStateHelper,
+                                                                   resv_coeff,
+                                                                   group_target->production_cmode};
     const auto& rates = ws.surface_rates;
     const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
     Scalar scale = 1.0;
     if (current_rate > 1e-14)
-        scale = target_rate->target_value / current_rate;
+        scale = group_target->target_value / current_rate;
 
     return scale;
 }

opm/simulators/wells/WellInterface.hpp

@@ -177,6 +177,10 @@ public:
                                         WellStateType& well_state,
                                         const bool solving_with_zero_rate);
 
+    // Convenience overload that gets scaled fractions internally
+    void updateGroupTargetFallbackFlag(WellStateType& well_state,
+                                       DeferredLogger& deferred_logger) const;
+
     // TODO: better name or further refactoring the function to make it more clear
     void prepareWellBeforeAssembling(const Simulator& simulator,
                                      const double dt,
@@ -320,6 +324,11 @@ public:
     computeCurrentWellRates(const Simulator& simulator,
                             DeferredLogger& deferred_logger) const = 0;
 
+    /// Get scaled well fractions from primary variables
+    /// Implemented by derived classes that have access to primary_variables_
+    virtual void getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
+                                        DeferredLogger& deferred_logger) const = 0;
+
     /// Modify the well_state's rates if there is only one nonzero rate.
     /// If so, that rate is kept as is, but the others are set proportionally
     /// to the rates at bhp limit or bhp 1 bar.

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -953,6 +953,72 @@ onlyKeepBHPandTHPcontrols(const SummaryState& summary_state,
     }
 }
 
+template<typename Scalar, typename IndexTraits>
+void WellInterfaceGeneric<Scalar, IndexTraits>::
+updateGroupTargetFallbackFlag(WellState<Scalar, IndexTraits>& well_state,
+                              const std::vector<Scalar>& scaled_well_fractions,
+                              DeferredLogger& deferred_logger) const
+{
+    auto& ws = well_state.well(this->index_of_well_);
+    if (!wellUnderGroupControl(ws) || this->isInjector() || 
+        !ws.group_target_fallback.has_value() || !ws.group_target.has_value())
+    {
+        ws.use_group_target_fallback = false;
+        return;
+    }
+    // We only need to check modes that can result in small fractions
+    // i.e. ORAT, WRAT, GRAT, LRAT
+    const auto cmode = ws.group_target->production_cmode;
+    const bool do_check = (cmode == Group::ProductionCMode::ORAT ||
+                           cmode == Group::ProductionCMode::WRAT ||
+                           cmode == Group::ProductionCMode::GRAT ||
+                           cmode == Group::ProductionCMode::LRAT);
+    if (!do_check) {
+        ws.use_group_target_fallback = false;
+        return;
+    }
+    // Check whether guiderate ratio or well fraction of production_cmode is too small, 
+    // if so, switch to fallback (original) control mode
+    const Scalar fraction_tolerance = this->param_.group_control_fraction_tolerance_;
+    if (ws.group_target->guiderate_ratio < fraction_tolerance) {
+        if (ws.group_target_fallback->guiderate_ratio < fraction_tolerance) {
+            // both current and fallback guiderate ratios are too small, don't switch
+            // This situation can be result of "badly" setup schedule.
+            const auto msg = fmt::format("Well {} guiderate ratio for current group control {} is {}, and fallback control {} is {},"
+                               " both are smaller than {}, we will keep current group control mode {}",
+                               name(), Group::ProductionCMode2String(cmode), ws.group_target->guiderate_ratio,
+                               Group::ProductionCMode2String(ws.group_target_fallback->production_cmode), ws.group_target_fallback->guiderate_ratio,
+                               fraction_tolerance, Group::ProductionCMode2String(cmode));
+            deferred_logger.info(msg);
+            ws.use_group_target_fallback = false;
+        } else {
+            ws.use_group_target_fallback = true;
+        }
+        return;
+    }
+    // check well fraction corresponding to cmode
+    Scalar cmode_frac = 0.0;
+    const auto& pu = this->phaseUsage();
+    if (cmode == Group::ProductionCMode::ORAT || cmode == Group::ProductionCMode::LRAT) {
+        if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) {
+            const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx);
+            cmode_frac += scaled_well_fractions[oil_pos];
+        }
+    }
+    if (cmode == Group::ProductionCMode::WRAT || cmode == Group::ProductionCMode::LRAT) {
+        if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) {
+            const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx);
+            cmode_frac += scaled_well_fractions[water_pos];
+        }
+    }
+    if (cmode == Group::ProductionCMode::GRAT) {
+        if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) {
+            const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx);
+            cmode_frac += scaled_well_fractions[gas_pos];
+        }
+    }
+    ws.use_group_target_fallback = (cmode_frac < fraction_tolerance); 
+}
 
 template class WellInterfaceGeneric<double, BlackOilDefaultFluidSystemIndices>;
 

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -249,6 +249,10 @@ protected:
         std::ranges::fill(this->inj_multiplier_damp_factor_, 1.0);
     }
 
+    void updateGroupTargetFallbackFlag(WellState<Scalar, IndexTraits>& well_state,
+                                       const std::vector<Scalar>& scaled_well_fractions,
+                                       DeferredLogger& deferred_logger) const;
+
     // definition of the struct OperabilityStatus
     struct OperabilityStatus
     {

opm/simulators/wells/WellInterface_impl.hpp

@@ -591,6 +591,15 @@ namespace Opm
                 }
             }
         }
+        // Add debug info for problematic group targets
+        if (this->isProducer() &&  ws.production_cmode == Well::ProducerCMode::GRUP && ws.use_group_target_fallback) {
+            assert(ws.group_target && ws.group_target_fallback);
+            const std::string cmode = Group::ProductionCMode2String(ws.group_target->production_cmode);
+            const std::string cmode_fallback = Group::ProductionCMode2String(ws.group_target_fallback->production_cmode);
+            const auto msg = fmt::format("    Well {} was solved using group target fallback mode {} as current group mode {} was not feasible.",
+                                            this->name(), cmode_fallback, cmode);
+            deferred_logger.debug(msg);
+        }
 
         return converged;
     }
@@ -988,6 +997,21 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateGroupTargetFallbackFlag(WellStateType& well_state,
+                                   DeferredLogger& deferred_logger) const
+    {
+        // Get scaled well fractions from derived class
+        std::vector<Scalar> scaled_well_fractions(FluidSystem::numPhases, 0.0);
+        this->getScaledWellFractions(scaled_well_fractions, deferred_logger);
+        // Call the base class method
+        this->Base::updateGroupTargetFallbackFlag(well_state, scaled_well_fractions, deferred_logger);
+    }
+
+
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -1664,6 +1688,7 @@ namespace Opm
             case Well::ProducerCMode::GRUP:
             {
                 assert(well.isAvailableForGroupControl());
+                this->updateGroupTargetFallbackFlag(well_state, deferred_logger);
                 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
                 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
                                                 well_state[well.name()].efficiency_scaling_factor;

regressionTests.cmake

@@ -813,22 +813,6 @@ add_test_compareECLFiles(CASENAME nonnc
                          REL_TOL ${rel_tol}
                          DIR editnnc)
 
-add_test_compareECLFiles(CASENAME nnc_overlapping_editnnc
-                         FILENAME NNC_OVERLAPPING_EDITNNC
-                         SIMULATOR flow
-                         ABS_TOL ${abs_tol}
-                         REL_TOL ${rel_tol}
-                         DIR editnnc
-                         TEST_ARGS --enable-opm-rst-file=1)
-
-add_test_compareECLFiles(CASENAME nnc_overlapping_editnncr_multregt
-                         FILENAME NNC_OVERLAPPING_EDITNNCR_MULTREGT
-                         SIMULATOR flow
-                         ABS_TOL ${abs_tol}
-                         REL_TOL ${rel_tol}
-                         DIR editnnc
-                         TEST_ARGS --enable-opm-rst-file=1)
-
 add_test_compareECLFiles(CASENAME spe1_foam
                          FILENAME SPE1FOAM
                          SIMULATOR flow