Merge pull request #4359 from totto82/fixGPMAINT2

Fix GPMAINT for groups without control

opm/simulators/wells/BlackoilWellModel.hpp

@@ -347,7 +347,7 @@ namespace Opm {
             bool alternative_well_rate_init_{};
 
             std::unique_ptr<RateConverterType> rateConverter_{};
-            std::unique_ptr<AverageRegionalPressureType> regionalAveragePressureCalculator_{};
+            std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};
 
 
             SimulatorReportSingle last_report_{};

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -208,9 +208,8 @@ namespace Opm {
 
             // Compute regional average pressures used by gpmaint
             if (schedule_[timeStepIdx].has_gpmaint()) {
-                const auto& fp = this->eclState_.fieldProps();
+                WellGroupHelpers::setRegionAveragePressureCalculator(fieldGroup, schedule(),
-                const auto& fipnum = fp.get_int("FIPNUM");
+                        timeStepIdx, this->eclState_.fieldProps(), phase_usage_, regionalAveragePressureCalculator_);
-                regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum));
             }
 
             {
@@ -330,10 +329,12 @@ namespace Opm {
         auto exc_type = ExceptionType::NONE;
         // update gpmaint targets
         if (schedule_[reportStepIdx].has_gpmaint()) {
-            regionalAveragePressureCalculator_->template defineState<ElementContext>(ebosSimulator_);
+            for (auto& calculator : regionalAveragePressureCalculator_) {
+                calculator.second->template defineState<ElementContext>(ebosSimulator_);
+            }
             const double dt = ebosSimulator_.timeStepSize();
             WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
-                                                           schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
+                                                           schedule_, regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
         }
         try {
             // Compute initial well solution for new wells and injectors that change injection type i.e. WAG.

opm/simulators/wells/TargetCalculator.cpp

@@ -95,29 +95,35 @@ RateType TargetCalculator::calcModeRateFromRates(const RateType* rates) const
     }
 }
 
-double TargetCalculator::groupTarget(const Group::ProductionControls ctrl) const
+double TargetCalculator::groupTarget(const std::optional<Group::ProductionControls>& ctrl, Opm::DeferredLogger& deferred_logger) const
 {
+    if (!ctrl && !use_gpmaint_) {
+        OPM_DEFLOG_THROW(std::logic_error,
+                         "Production group " + this->group_name_
+                         + "must either have a valid control or use GPMAINT",
+                         deferred_logger);
+    }
     switch (cmode_) {
     case Group::ProductionCMode::ORAT:
-        return ctrl.oil_target;
+        return ctrl->oil_target;
     case Group::ProductionCMode::WRAT:
-        return ctrl.water_target;
+        return ctrl->water_target;
     case Group::ProductionCMode::GRAT:
     {
         // gas target may have been adjusted by GCONSALE
         if ( group_grat_target_from_sales_ > 0)
             return group_grat_target_from_sales_;
 
-        return ctrl.gas_target;
+        return ctrl->gas_target;
     }
     case Group::ProductionCMode::LRAT:
-        return ctrl.liquid_target;
+        return ctrl->liquid_target;
     case Group::ProductionCMode::RESV:
     {
-        if(use_gpmaint_)
+        if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
             return this->group_state_.gpmaint_target(this->group_name_);
 
-        return ctrl.resv_target;
+        return ctrl->resv_target;
     }
     default:
         // Should never be here.
@@ -192,34 +198,40 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode
 }
 
 
-double InjectionTargetCalculator::groupTarget(const Group::InjectionControls& ctrl, Opm::DeferredLogger& deferred_logger) const
+double InjectionTargetCalculator::groupTarget(const std::optional<Group::InjectionControls>& ctrl, Opm::DeferredLogger& deferred_logger) const
 {
+    if (!ctrl && !use_gpmaint_) {
+        OPM_DEFLOG_THROW(std::logic_error,
+                         "Injection group " + this->group_name_
+                         + "must either have a valid control or use GPMAINT",
+                         deferred_logger);
+    }
     switch (cmode_) {
     case Group::InjectionCMode::RATE:
         if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
             return this->group_state_.gpmaint_target(this->group_name_);
 
-        return ctrl.surface_max_rate;
+        return ctrl->surface_max_rate;
     case Group::InjectionCMode::RESV:
         if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
             return this->group_state_.gpmaint_target(this->group_name_) / resv_coeff_[pos_];
 
-        return ctrl.resv_max_rate / resv_coeff_[pos_];
+        return ctrl->resv_max_rate / resv_coeff_[pos_];
     case Group::InjectionCMode::REIN: {
-        double production_rate = this->group_state_.injection_rein_rates(ctrl.reinj_group)[pos_];
+        double production_rate = this->group_state_.injection_rein_rates(ctrl->reinj_group)[pos_];
-        return ctrl.target_reinj_fraction * production_rate;
+        return ctrl->target_reinj_fraction * production_rate;
     }
     case Group::InjectionCMode::VREP: {
         const std::vector<double>& group_injection_reductions = this->group_state_.injection_reduction_rates(this->group_name_);
-        double voidage_rate = group_state_.injection_vrep_rate(ctrl.voidage_group) * ctrl.target_void_fraction;
+        double voidage_rate = group_state_.injection_vrep_rate(ctrl->voidage_group) * ctrl->target_void_fraction;
         double inj_reduction = 0.0;
-        if (ctrl.phase != Phase::WATER)
+        if (ctrl->phase != Phase::WATER)
             inj_reduction += group_injection_reductions[pu_.phase_pos[BlackoilPhases::Aqua]]
                     * resv_coeff_[pu_.phase_pos[BlackoilPhases::Aqua]];
-        if (ctrl.phase != Phase::OIL)
+        if (ctrl->phase != Phase::OIL)
             inj_reduction += group_injection_reductions[pu_.phase_pos[BlackoilPhases::Liquid]]
                     * resv_coeff_[pu_.phase_pos[BlackoilPhases::Liquid]];
-        if (ctrl.phase != Phase::GAS)
+        if (ctrl->phase != Phase::GAS)
             inj_reduction += group_injection_reductions[pu_.phase_pos[BlackoilPhases::Vapour]]
                     * resv_coeff_[pu_.phase_pos[BlackoilPhases::Vapour]];
         voidage_rate -= inj_reduction;

opm/simulators/wells/TargetCalculator.hpp

@@ -59,7 +59,7 @@ namespace WellGroupHelpers
         template <typename RateType>
         RateType calcModeRateFromRates(const RateType* rates) const;
 
-        double groupTarget(const Group::ProductionControls ctrl) const;
+        double groupTarget(const std::optional<Group::ProductionControls>& ctrl, Opm::DeferredLogger& deferred_logger) const;
 
         GuideRateModel::Target guideTargetMode() const;
 
@@ -94,7 +94,7 @@ namespace WellGroupHelpers
             return rates[pos_];
         }
 
-        double groupTarget(const Group::InjectionControls& ctrl, Opm::DeferredLogger& deferred_logger) const;
+        double groupTarget(const std::optional<Group::InjectionControls>& ctrl, Opm::DeferredLogger& deferred_logger) const;
 
         GuideRateModel::Target guideTargetMode() const;
 

opm/simulators/wells/WellAssemble.cpp

@@ -164,7 +164,8 @@ assembleControlEqProd(const WellState& well_state,
                                                              bhp, active_rates,
                                                              rCoeff,
                                                              efficiencyFactor,
-                                                             control_eq);
+                                                             control_eq,
+                                                             deferred_logger);
         break;
     }
     case Well::ProducerCMode::CMODE_UNDEFINED: {

opm/simulators/wells/WellGroupControls.cpp

@@ -155,8 +155,11 @@ getGroupInjectionControl(const Group& group,
         return tcalc.calcModeRateFromRates(groupTargetReductions);
     };
 
-    const double orig_target = tcalc.groupTarget(group.injectionControls(injectionPhase,
+    std::optional<Group::InjectionControls> ctrl;
-                                                                         summaryState),
+    if (!group.has_gpmaint_control(injectionPhase, currentGroupControl))
+        ctrl = group.injectionControls(injectionPhase, summaryState);
+
+    const double orig_target = tcalc.groupTarget(ctrl,
                                                 deferred_logger);
     const auto chain = WellGroupHelpers::groupChainTopBot(well_.name(), group.name(),
                                                           schedule, well_.currentStep());
@@ -214,7 +217,6 @@ getGroupInjectionTargetRate(const Group& group,
         // Should not be here.
         assert(false);
     }
-
     auto currentGroupControl = group_state.injection_control(group.name(), injectionPhase);
     if (currentGroupControl == Group::InjectionCMode::FLD ||
         currentGroupControl == Group::InjectionCMode::NONE) {
@@ -274,7 +276,12 @@ getGroupInjectionTargetRate(const Group& group,
         return tcalc.calcModeRateFromRates(groupTargetReductions);
     };
 
-    const double orig_target = tcalc.groupTarget(group.injectionControls(injectionPhase, summaryState), deferred_logger);
+    std::optional<Group::InjectionControls> ctrl;
+    if (!group.has_gpmaint_control(injectionPhase, currentGroupControl))
+        ctrl = group.injectionControls(injectionPhase, summaryState);
+
+    const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
+
     const auto chain = WellGroupHelpers::groupChainTopBot(well_.name(), group.name(), schedule, well_.currentStep());
     // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
     const size_t num_ancestors = chain.size() - 1;
@@ -301,7 +308,8 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
                                                   const std::vector<EvalWell>& rates,
                                                   const RateConvFunc& rateConverter,
                                                   double efficiencyFactor,
-                                                  EvalWell& control_eq) const
+                                                  EvalWell& control_eq,
+                                                  DeferredLogger& deferred_logger) const
 {
     const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
     if (currentGroupControl == Group::ProductionCMode::FLD ||
@@ -325,7 +333,7 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
             getGroupProductionControl(parent, well_state, group_state,
                                       schedule, summaryState, bhp,
                                       rates, rateConverter,
-                                      efficiencyFactor, control_eq);
+                                      efficiencyFactor, control_eq, deferred_logger);
             return;
         }
     }
@@ -372,7 +380,11 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
         return tcalc.calcModeRateFromRates(groupTargetReductions);
     };
 
-    const double orig_target = tcalc.groupTarget(group.productionControls(summaryState));
+    std::optional<Group::ProductionControls> ctrl;
+    if (!group.has_gpmaint_control(currentGroupControl))
+        ctrl = group.productionControls(summaryState);
+
+    const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
     const auto chain = WellGroupHelpers::groupChainTopBot(well_.name(), group.name(),
                                                           schedule, well_.currentStep());
     // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
@@ -400,7 +412,8 @@ getGroupProductionTargetRate(const Group& group,
                              const Schedule& schedule,
                              const SummaryState& summaryState,
                              const RateConvFunc& rateConverter,
-                             double efficiencyFactor) const
+                             double efficiencyFactor,
+                             DeferredLogger& deferred_logger) const
 {
     const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
     if (currentGroupControl == Group::ProductionCMode::FLD ||
@@ -413,7 +426,8 @@ getGroupProductionTargetRate(const Group& group,
             efficiencyFactor *= group.getGroupEfficiencyFactor();
             return getGroupProductionTargetRate(parent, well_state, group_state,
                                                 schedule, summaryState,
-                                                rateConverter, efficiencyFactor);
+                                                rateConverter, efficiencyFactor,
+                                                deferred_logger);
         }
     }
 
@@ -452,7 +466,11 @@ getGroupProductionTargetRate(const Group& group,
         return tcalc.calcModeRateFromRates(groupTargetReductions);
     };
 
-    const double orig_target = tcalc.groupTarget(group.productionControls(summaryState));
+    std::optional<Group::ProductionControls> ctrl;
+    if (!group.has_gpmaint_control(currentGroupControl))
+        ctrl = group.productionControls(summaryState);
+
+    const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
     const auto chain = WellGroupHelpers::groupChainTopBot(well_.name(), group.name(),
                                                           schedule, well_.currentStep());
     // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
@@ -506,7 +524,8 @@ getGroupProductionControl<__VA_ARGS__>(const Group&, \
                                        const std::vector<__VA_ARGS__>&, \
                                        const RateConvFunc& rateConverter, \
                                        double efficiencyFactor, \
-                                       __VA_ARGS__& control_eq) const;
+                                       __VA_ARGS__& control_eq, \
+                                       DeferredLogger& deferred_logger) const; \
 
 INSTANCE(DenseAd::Evaluation<double,3,0u>)
 INSTANCE(DenseAd::Evaluation<double,4,0u>)

opm/simulators/wells/WellGroupControls.hpp

@@ -84,7 +84,8 @@ public:
                                    const std::vector<EvalWell>& rates,
                                    const RateConvFunc& rateConverter,
                                    double efficiencyFactor,
-                                   EvalWell& control_eq) const;
+                                   EvalWell& control_eq,
+                                   DeferredLogger& deferred_logger) const;
 
     double getGroupProductionTargetRate(const Group& group,
                                         const WellState& well_state,
@@ -92,7 +93,8 @@ public:
                                         const Schedule& schedule,
                                         const SummaryState& summaryState,
                                         const RateConvFunc& rateConverter,
-                                        double efficiencyFactor) const;
+                                        double efficiencyFactor,
+                                        DeferredLogger& deferred_logger) const;
 
 private:
     const WellInterfaceGeneric& well_; //!< Reference to well interface

opm/simulators/wells/WellGroupHelpers.cpp

@@ -689,7 +689,7 @@ namespace WellGroupHelpers
             return;
 
         const auto [name, number] = *region;
-        const double error = gpm->pressure_target() - regional_values.pressure(number);
+        const double error = gpm->pressure_target() - regional_values.at(name)->pressure(number);
         double current_rate = 0.0;
         const auto& pu = well_state.phaseUsage();
         double sign = 1.0;
@@ -744,7 +744,6 @@ namespace WellGroupHelpers
             default:
                 throw std::invalid_argument("Invalid Flow target type in GPMAINT");
         }
-
         auto& gpmaint_state = group_state.gpmaint(group.name());
         double rate = gpm->rate(gpmaint_state, current_rate, error, dt);
         group_state.update_gpmaint_target(group.name(), std::max(0.0, sign * rate));
@@ -1268,7 +1267,12 @@ namespace WellGroupHelpers
             const std::vector<double>& groupSurfaceRates = group_state.production_rates(group_name);
             return tcalc.calcModeRateFromRates(groupSurfaceRates);
         };
-        const double orig_target = tcalc.groupTarget(group.productionControls(summaryState));
+
+        std::optional<Group::ProductionControls> ctrl;
+        if (!group.has_gpmaint_control(currentGroupControl))
+            ctrl = group.productionControls(summaryState);
+
+        const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
         // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP.
         // Then ...
         // TODO finish explanation.
@@ -1410,7 +1414,12 @@ namespace WellGroupHelpers
             return tcalc.calcModeRateFromRates(groupSurfaceRates);
         };
 
-        const double orig_target = tcalc.groupTarget(group.injectionControls(injectionPhase, summaryState), deferred_logger);
+        std::optional<Group::InjectionControls> ctrl;
+        if (!group.has_gpmaint_control(injectionPhase, currentGroupControl))
+            ctrl = group.injectionControls(injectionPhase, summaryState);
+
+
+        const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
         // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP.
         // Then ...
         // TODO finish explanation.
@@ -1477,6 +1486,32 @@ namespace WellGroupHelpers
         return std::make_pair(current_rate > target_rate, scale);
     }
 
+    template <class AverageRegionalPressureType>
+    void setRegionAveragePressureCalculator(const Group& group,
+                                            const Schedule& schedule,
+                                            const int reportStepIdx,
+                                            const FieldPropsManager& fp,
+                                            const PhaseUsage& pu,
+                                            std::map<std::string, std::unique_ptr<AverageRegionalPressureType>>& regionalAveragePressureCalculator) 
+    {
+        for (const std::string& groupName : group.groups()) {
+            setRegionAveragePressureCalculator( schedule.getGroup(groupName, reportStepIdx), schedule, 
+                                                reportStepIdx, fp, pu, regionalAveragePressureCalculator);
+        }
+        const auto& gpm = group.gpmaint();
+        if (!gpm)
+            return;
+
+        const auto& reg = gpm->region();
+        if (!reg)
+            return;
+
+        if (regionalAveragePressureCalculator.count(reg->first) == 0) {
+            const std::string name = (reg->first.rfind("FIP", 0) == 0) ? reg->first : "FIP" + reg->first;
+            const auto& fipnum = fp.get_int(name);
+            regionalAveragePressureCalculator[reg->first] = std::make_unique<AverageRegionalPressureType>(pu,fipnum);
+        }
+    }
 
     template <class Comm>
     void updateGuideRates(const Group& group,
@@ -1615,13 +1650,20 @@ namespace WellGroupHelpers
     }
 
  using AvgP = RegionAverageCalculator::AverageRegionalPressure<BlackOilFluidSystem<double>,std::vector<int>>;
- template void WellGroupHelpers::updateGpMaintTargetForGroups<AvgP>(const Group&,
+ using AvgPMap = std::map<std::string, std::unique_ptr<AvgP>>;
+ template void WellGroupHelpers::updateGpMaintTargetForGroups<AvgPMap>(const Group&,
                                                                     const Schedule&,
-                                                                    const AvgP&,
+                                                                    const AvgPMap&,
                                                                     int,
                                                                     double,
                                                                     const WellState&,
                                                                     GroupState&);
+ template void WellGroupHelpers::setRegionAveragePressureCalculator<AvgP>(const Group&,
+                                            const Schedule&,
+                                            const int,
+                                            const FieldPropsManager&,
+                                            const PhaseUsage&,
+                                            AvgPMap&);
 
 template
 void updateGuideRateForProductionGroups<Parallel::Communication>(const Group& group,

opm/simulators/wells/WellGroupHelpers.hpp

@@ -22,6 +22,8 @@
 #define OPM_WELLGROUPHELPERS_HEADER_INCLUDED
 
 #include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
+
 
 #include <map>
 #include <string>
@@ -38,6 +40,7 @@ struct PhaseUsage;
 class Schedule;
 class VFPProdProperties;
 class WellState;
+class FieldPropsManager;
 
 namespace Network { class ExtNetwork; }
 
@@ -309,6 +312,14 @@ namespace WellGroupHelpers
                                                       const std::vector<double>& resv_coeff,
                                                       DeferredLogger& deferred_logger);
 
+    template <class AverageRegionalPressureType>
+    void setRegionAveragePressureCalculator(const Group& group,
+                                            const Schedule& schedule,
+                                            const int reportStepIdx,
+                                            const FieldPropsManager& fp,
+                                            const PhaseUsage& pu,
+                                            std::map<std::string, std::unique_ptr<AverageRegionalPressureType>>& regionalAveragePressureCalculator);
+
 
 } // namespace WellGroupHelpers
 

opm/simulators/wells/WellInterfaceFluidSystem.cpp

@@ -189,7 +189,8 @@ getGroupProductionTargetRate(const Group& group,
                           const GroupState& group_state,
                           const Schedule& schedule,
                           const SummaryState& summaryState,
-                          double efficiencyFactor) const
+                          double efficiencyFactor,
+                          DeferredLogger& deferred_logger) const
 {
     auto rCoeff = [this](const int id, const int region, std::vector<double>& coeff)
     {
@@ -199,7 +200,8 @@ getGroupProductionTargetRate(const Group& group,
     return WellGroupControls(*this).getGroupProductionTargetRate(group, well_state,
                                                                  group_state, schedule,
                                                                  summaryState,
-                                                                 rCoeff, efficiencyFactor);
+                                                                 rCoeff, efficiencyFactor,
+                                                                 deferred_logger);
 }
 
 template class WellInterfaceFluidSystem<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;

opm/simulators/wells/WellInterfaceFluidSystem.hpp

@@ -109,7 +109,8 @@ protected:
                                  const GroupState& group_state,
                                  const Schedule& schedule,
                                  const SummaryState& summaryState,
-                                 double efficiencyFactor) const;
+                                 double efficiencyFactor,
+                                 DeferredLogger& deferred_logger) const;
 
     // For the conversion between the surface volume rate and reservoir voidage rate
     const RateConverterType& rateConverter_;

opm/simulators/wells/WellInterface_impl.hpp

@@ -997,7 +997,8 @@ namespace Opm
                                                           group_state,
                                                           schedule,
                                                           summaryState,
-                                                          efficiencyFactor);
+                                                          efficiencyFactor,
+                                                          deferred_logger);
 
                 // we don't want to scale with zero and get zero rates.
                 if (scale > 0) {