Merge pull request #6837 from hakonhagland/refact_wgc

Move getAutoChokeGroupProductionTargetRate() to GroupStateHelper

opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp

@@ -38,7 +38,6 @@
 #include <opm/simulators/wells/BlackoilWellModel.hpp>
 #include <opm/simulators/wells/TargetCalculator.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
-#include <opm/simulators/wells/WellGroupControls.hpp>
 
 #include <fmt/format.h>
 
@@ -192,10 +191,9 @@ computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger)
                 // derived via group guide rates
                 const Scalar efficiencyFactor = 1.0;
                 const Group& parentGroup = well_model_.schedule().getGroup(group.parent(), reportStepIdx);
-                auto target = WellGroupControls<Scalar, IndexTraits>::
+                auto target = well_model_.groupStateHelper().
                     getAutoChokeGroupProductionTargetRate(group,
                                                           parentGroup,
-                                                          well_model_.groupStateHelper(),
                                                           resv_coeff,
                                                           efficiencyFactor);
                 target_tmp = target.first;

opm/simulators/wells/GroupStateHelper.cpp

@@ -189,7 +189,7 @@ GroupStateHelper<Scalar, IndexTraits>::checkGroupConstraintsInj(const std::strin
 
     const Scalar orig_target = tcalc.groupTarget();
     const Scalar current_rate_available = tcalc.calcModeRateFromRates(rates);
-    const std::size_t local_reduction_level = this->getLocalReductionLevel(
+    const std::size_t local_reduction_level = this->getLocalReductionLevel_(
         chain, /*is_production_group=*/false, injection_phase);
 
     const Scalar target = this->applyReductionsAndFractions_(chain,
@@ -330,7 +330,7 @@ GroupStateHelper<Scalar, IndexTraits>::checkGroupConstraintsProd(const std::stri
         return std::make_pair(current_well_rate_available > group_target_rate_available, scale);
     }
 
-    const std::size_t local_reduction_level = this->getLocalReductionLevel(
+    const std::size_t local_reduction_level = this->getLocalReductionLevel_(
         chain, /*is_production_group=*/true, /*injection_phase=*/Phase::OIL);
 
     const Scalar target = this->applyReductionsAndFractions_(chain,
@@ -515,6 +515,82 @@ getInjectionGuideTargetMode(Phase injection_phase) const
     }
 }
 
+template<typename Scalar, typename IndexTraits>
+std::pair<Scalar, Group::ProductionCMode>
+GroupStateHelper<Scalar, IndexTraits>::
+getAutoChokeGroupProductionTargetRate(const Group& bottom_group,
+                                      const Group& group,
+                                      const std::vector<Scalar>& resv_coeff,
+                                      Scalar efficiencyFactor) const
+{
+    const Group::ProductionCMode& currentGroupControl
+        = this->groupState().production_control(group.name());
+    if (currentGroupControl == Group::ProductionCMode::FLD ||
+        currentGroupControl == Group::ProductionCMode::NONE) {
+        if (!group.productionGroupControlAvailable()) {
+            return std::make_pair(1.0, currentGroupControl);
+        } else {
+            // Produce share of parents control
+            const auto& parent = this->schedule_.getGroup(group.parent(), this->report_step_);
+            efficiencyFactor *= group.getGroupEfficiencyFactor();
+            return this->getAutoChokeGroupProductionTargetRate(bottom_group, parent,
+                                                               resv_coeff, efficiencyFactor);
+        }
+    }
+
+    if (!group.isProductionGroup()) {
+        return std::make_pair(1.0, currentGroupControl);
+    }
+
+    // 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.
+
+    GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{*this,
+                                                                    resv_coeff,
+                                                                    group};
+
+    GroupStateHelpers::FractionCalculator<Scalar, IndexTraits> fcalc(this->schedule_,
+                                                                      *this,
+                                                                      this->summary_state_,
+                                                                      this->report_step_,
+                                                                      &this->guide_rate_,
+                                                                      tcalc.guideTargetMode(),
+                                                                      true,
+                                                                      Phase::OIL);
+
+    auto localFraction = [&](const std::string& child) {
+        return fcalc.localFraction(child, child); //Note child needs to be passed to always include since the global isGrup map is not updated yet.
+    };
+
+    auto localReduction = [&](const std::string& group_name) {
+        const std::vector<Scalar>& groupTargetReductions
+            = this->groupState().production_reduction_rates(group_name);
+        return tcalc.calcModeRateFromRates(groupTargetReductions);
+    };
+
+    const Scalar orig_target = tcalc.groupTarget();
+    const auto chain = this->groupChainTopBot(bottom_group.name(), group.name());
+    const std::size_t local_reduction_level = this->getLocalReductionLevel_(
+        chain, /*is_production_group=*/true, Phase::OIL);
+
+    // Delegate to applyReductionsAndFractions_ with no addback.  Addback is not
+    // needed because autochoke wells are excluded from reductions via the
+    // !group.as_choke() check in updateGroupTargetReductionRecursive_().
+    const Scalar target = this->applyReductionsAndFractions_(chain,
+                                                              orig_target,
+                                                              /*current_rate_available=*/Scalar(0.0),
+                                                              local_reduction_level,
+                                                              /*is_production_group=*/true,
+                                                              /*injection_phase=*/Phase::OIL,
+                                                              localReduction,
+                                                              localFraction,
+                                                              /*do_addback=*/false);
+    // Avoid negative target rates coming from too large local reductions.
+    const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor);
+
+    return std::make_pair(target_rate, currentGroupControl);
+}
+
 template <typename Scalar, typename IndexTraits>
 std::vector<Scalar>
 GroupStateHelper<Scalar, IndexTraits>::
@@ -615,7 +691,7 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::str
     const Scalar current_rate_available = tcalc.calcModeRateFromRates(rates);
     const auto chain = this->groupChainTopBot(name, group.name());
 
-    const std::size_t local_reduction_level = this->getLocalReductionLevel(
+    const std::size_t local_reduction_level = this->getLocalReductionLevel_(
         chain, /*is_production_group=*/false, injection_phase);
 
     // Add-back is only performed when the well is under individual control (not GRUP)
@@ -711,7 +787,7 @@ GroupStateHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::str
     const Scalar current_rate_available = -tcalc.calcModeRateFromRates(rates);
     const auto chain = this->groupChainTopBot(name, group.name());
 
-    const std::size_t local_reduction_level = this->getLocalReductionLevel(
+    const std::size_t local_reduction_level = this->getLocalReductionLevel_(
         chain, /*is_production_group=*/true, /*injection_phase=*/Phase::OIL);
 
     // Add-back is only performed when the well is under individual control (not GRUP)
@@ -1540,7 +1616,7 @@ GroupStateHelper<Scalar, IndexTraits>::getGuideRateVector_(const std::vector<Sca
 
 template <typename Scalar, typename IndexTraits>
 std::size_t
-GroupStateHelper<Scalar, IndexTraits>::getLocalReductionLevel(const std::vector<std::string>& chain,
+GroupStateHelper<Scalar, IndexTraits>::getLocalReductionLevel_(const std::vector<std::string>& chain,
                                                                const bool is_production_group,
                                                                const Phase injection_phase) const
 {
@@ -1630,14 +1706,14 @@ GroupStateHelper<Scalar, IndexTraits>::isAutoChokeGroupUnderperforming_(const Gr
 
     // Calculate target rate for this auto choke group using a guide rate ratio
     // formula. This is an approximation that ignores intermediate reductions.
-    // A more accurate calculation would use WellGroupControls::getAutoChokeGroupProductionTargetRate(),
+    // A more accurate calculation would use getAutoChokeGroupProductionTargetRate(),
     // but that requires reduction rates not yet available during
     // updateGroupControlledWells() (circular dependency). The efficiency factor is however accounted for
     // below meaning that the current target rate will be correct for the special case where the control
     // group has no reductions and the auto choke group is a direct child of the control group.
     // TODO: The issue with missing handling of reductions should be fixed by refactoring the auto choke
     // target calculation below to use the more accurate
-    // WellGroupControls::getAutoChokeGroupProductionTargetRate() function.
+    // getAutoChokeGroupProductionTargetRate() function.
     const auto num_phases = this->numPhases();
     std::vector<Scalar> resv_coeff(num_phases, 1.0);
     GroupStateHelpers::TargetCalculator<Scalar, IndexTraits> tcalc{*this,

opm/simulators/wells/GroupStateHelper.hpp

@@ -243,20 +243,6 @@ public:
     /// @return The corresponding GuideRateModel::Target for the injection phase
     GuideRateModel::Target getInjectionGuideTargetMode(Phase injection_phase) const;
 
-    //! \brief Find the local reduction level in a group chain.
-    //!
-    //! The local reduction level is the deepest level in the chain (starting from level 1)
-    //! where a group has both a guide rate and group-controlled wells (GCW > 0).
-    //! This level determines where reductions are applied and add-back occurs.
-    //!
-    //! \param chain Group chain from control group (top) to bottom group
-    //! \param is_production_group True for production, false for injection
-    //! \param injection_phase Phase for injection groups (ignored for production)
-    //! \return The local reduction level (0 if no intermediate group qualifies)
-    std::size_t getLocalReductionLevel(const std::vector<std::string>& chain,
-        bool is_production_group,
-        Phase injection_phase) const;
-
     Scalar getProductionGroupTarget(const Group& group) const;
 
     /// @brief Get the guide rate target mode for a production group
@@ -264,6 +250,12 @@ public:
     /// @return The GuideRateModel::Target based on the group's production control mode
     GuideRateModel::Target getProductionGuideTargetMode(const Group& group) const;
 
+    std::pair<Scalar, Group::ProductionCMode>
+    getAutoChokeGroupProductionTargetRate(const Group& bottom_group,
+                                          const Group& group,
+                                          const std::vector<Scalar>& resv_coeff,
+                                          Scalar efficiencyFactor) const;
+
     GuideRate::RateVector getProductionGroupRateVector(const std::string& group_name) const;
 
     using GroupTarget = typename SingleWellState<Scalar, IndexTraits>::GroupTarget;
@@ -537,6 +529,19 @@ private:
 
     GuideRate::RateVector getGuideRateVector_(const std::vector<Scalar>& rates) const;
 
+    //! \brief Find the local reduction level in a group chain.
+    //!
+    //! The local reduction level is the deepest level in the chain (starting from level 1)
+    //! where a group has both a guide rate and group-controlled wells (GCW > 0).
+    //! This level determines where reductions are applied and add-back occurs.
+    //!
+    //! \param chain Group chain from control group (top) to bottom group
+    //! \param is_production_group True for production, false for injection
+    //! \param injection_phase Phase for injection groups (ignored for production)
+    //! \return The local reduction level (0 if no intermediate group qualifies)
+    std::size_t getLocalReductionLevel_(const std::vector<std::string>& chain,
+        bool is_production_group,
+        Phase injection_phase) const;
 
 #ifdef RESERVOIR_COUPLING_ENABLED
     Scalar getReservoirCouplingMasterGroupRate_(const Group& group,

opm/simulators/wells/WellGroupControls.cpp

@@ -31,7 +31,6 @@
 #include <opm/material/densead/Evaluation.hpp>
 #include <opm/material/fluidsystems/BlackOilDefaultFluidSystemIndices.hpp>
 
-#include <opm/simulators/wells/FractionCalculator.hpp>
 #include <opm/simulators/wells/GroupState.hpp>
 #include <opm/simulators/wells/TargetCalculator.hpp>
 #include <opm/simulators/wells/GroupStateHelper.hpp>
@@ -328,92 +327,6 @@ getGroupProductionTargetRate(const Group& group,
     return scale;
 }
 
-template<typename Scalar, typename IndexTraits>
-std::pair<Scalar, Group::ProductionCMode> WellGroupControls<Scalar, IndexTraits>::
-getAutoChokeGroupProductionTargetRate(const Group& bottom_group,
-                                      const Group& group,
-                                      const GroupStateHelperType& groupStateHelper,
-                                      const std::vector<Scalar>& resv_coeff,
-                                      Scalar efficiencyFactor)
-{
-    const auto& schedule = groupStateHelper.schedule();
-    const auto& summaryState = groupStateHelper.summaryState();
-    const int reportStepIdx = groupStateHelper.reportStepIdx();
-    const GuideRate* guideRate = &groupStateHelper.guideRate();
-    const auto& group_state = groupStateHelper.groupState();
-    const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
-    if (currentGroupControl == Group::ProductionCMode::FLD ||
-        currentGroupControl == Group::ProductionCMode::NONE) {
-        if (!group.productionGroupControlAvailable()) {
-            return std::make_pair(1.0, currentGroupControl);
-        } else {
-            // Produce share of parents control
-            const auto& parent = schedule.getGroup(group.parent(), reportStepIdx);
-            efficiencyFactor *= group.getGroupEfficiencyFactor();
-            return getAutoChokeGroupProductionTargetRate(bottom_group, parent, groupStateHelper,
-                                                resv_coeff, efficiencyFactor);
-        }
-    }
-
-    if (!group.isProductionGroup()) {
-        return std::make_pair(1.0, currentGroupControl);
-    }
-
-    // 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.
-
-    // Make conversion factors for RESV <-> surface rates.
-    // std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 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};
-
-    GroupStateHelpers::FractionCalculator<Scalar, IndexTraits> fcalc(schedule,
-                                                                     groupStateHelper,
-                                                                     summaryState,
-                                                                     reportStepIdx,
-                                                                     guideRate,
-                                                                     tcalc.guideTargetMode(),
-                                                                     true,
-                                                                     Phase::OIL);
-
-    auto localFraction = [&](const std::string& child) {
-        return fcalc.localFraction(child, child); //Note child needs to be passed to always include since the global isGrup map is not updated yet.
-    };
-
-    auto localReduction = [&](const std::string& group_name) {
-        const std::vector<Scalar>& groupTargetReductions = group_state.production_reduction_rates(group_name);
-        return tcalc.calcModeRateFromRates(groupTargetReductions);
-    };
-
-    const Scalar orig_target = tcalc.groupTarget();
-    const auto chain = groupStateHelper.groupChainTopBot(bottom_group.name(), group.name());
-    // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
-    const std::size_t num_ancestors = chain.size() - 1;
-    const std::size_t local_reduction_level = groupStateHelper.getLocalReductionLevel(
-        chain, /*is_production_group=*/true, Phase::OIL);
-    Scalar target = orig_target;
-    for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
-        if ((ii == 0) || guideRate->has(chain[ii])) {
-            // Apply local reductions only at the control level (top) and for
-            // intermediate levels where we have both a guide rate and
-            // group-controlled wells (local_reduction_level >= ii).  Without
-            // this gating, reductions from a group with GCW=0 would be
-            // double-counted.
-            if (local_reduction_level >= ii) {
-                target -= localReduction(chain[ii]);
-            }
-        }
-        target *= localFraction(chain[ii+1]);
-    }
-    // Avoid negative target rates coming from too large local reductions.
-    const double target_rate = std::max(Scalar(0.0), target / efficiencyFactor);
-
-    return std::make_pair(target_rate, currentGroupControl);
-}
-
 #define INSTANTIATE(T,...)                                               \
     template void WellGroupControls<T, BlackOilDefaultFluidSystemIndices>::                                 \
         getGroupInjectionControl(const Group&,                           \

opm/simulators/wells/WellGroupControls.hpp

@@ -86,12 +86,6 @@ public:
                                         const RateConvFunc& rateConverter,
                                         Scalar efficiencyFactor) const;
 
-    static std::pair<Scalar, Group::ProductionCMode> getAutoChokeGroupProductionTargetRate(const Group& bottom_group,
-                                                        const Group& group,
-                                                        const GroupStateHelperType& groupStateHelper,
-                                                        const std::vector<Scalar>& resv_coeff,
-                                                        Scalar efficiencyFactor);
-
 private:
     const WellInterfaceGeneric<Scalar, IndexTraits>& well_; //!< Reference to well interface
 };

tests/test_group_higher_constraints.cpp

@@ -62,7 +62,6 @@
 #include <opm/simulators/wells/GroupState.hpp>
 #include <opm/simulators/wells/FractionCalculator.hpp>
 #include <opm/simulators/wells/TargetCalculator.hpp>
-#include <opm/simulators/wells/WellGroupControls.hpp>
 
 #include <boost/test/unit_test.hpp>
 
@@ -1057,16 +1056,14 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupProductionTargetRate)
 
     // Call getAutoChokeGroupProductionTargetRate()
     // This is how it's called in BlackoilWellModelNetwork:
-    //   getAutoChokeGroupProductionTargetRate(bottom_group=MANI, group=MANI, gsh, resv_coeff, 1.0)
+    //   gsh.getAutoChokeGroupProductionTargetRate(bottom_group=MANI, group=MANI, resv_coeff, 1.0)
     // The function recurses upward from MANI (FLD) → PLAT (FLD) → FIELD (ORAT=controlling).
     // During recursion: efficiencyFactor accumulates E_MANI * E_PLAT = 0.8 * 0.9 = 0.72.
     auto resv_coeff = resvCoeff();
 
-    using WGC = Opm::WellGroupControls<double, IndexTraits>;
-    auto [target_rate_si, control_mode] = WGC::getAutoChokeGroupProductionTargetRate(
+    auto [target_rate_si, control_mode] = gsh.getAutoChokeGroupProductionTargetRate(
         maniGroup(),      // bottom_group
         maniGroup(),      // group (starts recursion here, recurses up to controlling group)
-        gsh,
         resv_coeff,
         /*efficiencyFactor=*/1.0
     );
@@ -1218,11 +1215,9 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateWithUnderperformance)
 
     // Call getAutoChokeGroupProductionTargetRate
     auto resv_coeff = resvCoeff();
-    using WGC = Opm::WellGroupControls<double, IndexTraits>;
-    auto [target_rate_si, control_mode] = WGC::getAutoChokeGroupProductionTargetRate(
+    auto [target_rate_si, control_mode] = gsh.getAutoChokeGroupProductionTargetRate(
         maniGroup(),      // bottom_group
         maniGroup(),      // group (starts recursion here)
-        gsh,
         resv_coeff,
         /*efficiencyFactor=*/1.0
     );