implement proper guiderate support

Tested on simple problems including GUIDERAT

opm/simulators/wells/BlackoilWellModel.hpp

@@ -35,6 +35,7 @@
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/WellCollection.hpp>
@@ -272,6 +273,7 @@ namespace Opm {
             SimulatorReport last_report_;
 
             WellTestState wellTestState_;
+            std::unique_ptr<GuideRate> guideRate_;
 
             // used to better efficiency of calcuation
             mutable BVector scaleAddRes_;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -34,6 +34,9 @@ namespace Opm {
         terminal_output_ = false;
         if (ebosSimulator.gridView().comm().rank() == 0)
             terminal_output_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
+
+        // Create the guide rate container.
+        guideRate_.reset(new GuideRate (ebosSimulator_.vanguard().schedule()));
     }
 
     template<typename TypeTag>
@@ -377,12 +380,33 @@ namespace Opm {
 
         for (auto& well : well_container_) {
             well->setVFPProperties(vfp_properties_.get());
+            well->setGuideRate(guideRate_.get());
         }
 
         // Close completions due to economical reasons
         for (auto& well : well_container_) {
             well->closeCompletions(wellTestState_);
         }
+
+        // calculate the well potentials
+        try {
+            std::vector<double> well_potentials;
+            computeWellPotentials(well_potentials, reportStepIdx, local_deferredLogger);
+        } catch ( std::runtime_error& e ) {
+            const std::string msg = "A zero well potential is returned for output purposes. ";
+            local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
+        }
+
+        //compute well guideRates
+        const int np = numPhases();
+        for (const auto& well : well_container_) {
+            const double& oilpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Liquid]];
+            const double& gaspot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Vapour]];
+            const double& waterpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Aqua]];
+            guideRate_->compute(well->name(), reportStepIdx, simulationTime, oilpot, gaspot, waterpot);
+        }
+        const Group2& fieldGroup = schedule().getGroup2("FIELD", reportStepIdx);
+        wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, guideRate_.get(), well_state_);
     }
 
 
@@ -412,6 +436,7 @@ namespace Opm {
                 wellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup2(wellEcl.groupName(), timeStepIdx), schedule(), timeStepIdx, well_efficiency_factor);
                 well->setWellEfficiencyFactor(well_efficiency_factor);
                 well->setVFPProperties(vfp_properties_.get());
+                well->setGuideRate(guideRate_.get());
 
                 const WellTestConfig::Reason testing_reason = testWell.second;
 
@@ -451,10 +476,10 @@ namespace Opm {
         // update the rate converter with current averages pressures etc in
         rateConverter_->template defineState<ElementContext>(ebosSimulator_);
 
+        const int reportStepIdx = ebosSimulator_.episodeIndex();
         // calculate the well potentials
         try {
             std::vector<double> well_potentials;
-            const int reportStepIdx = ebosSimulator_.episodeIndex();
             computeWellPotentials(well_potentials, reportStepIdx, local_deferredLogger);
         } catch ( std::runtime_error& e ) {
             const std::string msg = "A zero well potential is returned for output purposes. ";
@@ -467,6 +492,8 @@ namespace Opm {
             global_deferredLogger.logMessages();
         }
 
+
+
     }
 
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1811,31 +1811,39 @@ namespace Opm
         const auto& groupcontrols = group.injectionControls(summaryState);
 
         int phasePos;
+        Well2::GuideRateTarget wellTarget;
+        Group2::GuideRateTarget groupTarget;
+
         switch (injectorType) {
         case Well2::InjectorType::WATER:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Aqua];
+            wellTarget = Well2::GuideRateTarget::WAT;
+            groupTarget = Group2::GuideRateTarget::WAT;
             break;
         }
         case Well2::InjectorType::OIL:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Liquid];
+            wellTarget = Well2::GuideRateTarget::OIL;
+            groupTarget = Group2::GuideRateTarget::OIL;
             break;
         }
         case Well2::InjectorType::GAS:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Vapour];
+            wellTarget = Well2::GuideRateTarget::GAS;
+            groupTarget = Group2::GuideRateTarget::GAS;
             break;
         }
         default:
             throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
         }
 
-        double groupTotalGuideRate = 0.0;
         double groupTargetReduction = 0.0;
-        wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTotalGuideRate, groupTargetReduction);
-        double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-        double fraction = guideRate / groupTotalGuideRate;
+        wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTargetReduction);
+        double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, wellTarget, /*isInjector*/true);
+        wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, groupTarget, /*isInjector*/true, fraction);
 
         switch(currentGroupControl) {
         case Group2::InjectionCMode::NONE:
@@ -1924,12 +1932,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::ORAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::OIL, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::OIL, /*isInjector*/false, fraction);
 
             const double rate_target = std::max(0.0, groupcontrols.oil_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
@@ -1939,12 +1946,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::WRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::WAT, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::WAT, /*isInjector*/false, fraction);
 
             const double rate_target = std::max(0.0, groupcontrols.water_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
@@ -1954,13 +1960,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::GRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
-
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::GAS, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::GAS, /*isInjector*/false, fraction);
             const double rate_target = std::max(0.0, groupcontrols.gas_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::gasCompIdx));
             const EvalWell& rate = -getSegmentRate(0, Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
@@ -1969,15 +1973,14 @@ namespace Opm
         }
         case Group2::ProductionCMode::LRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
             phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            guideRate += wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::LIQ, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::LIQ, /*isInjector*/false, fraction);
+
             const double rate_target = std::max(0.0, groupcontrols.liquid_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -966,31 +966,40 @@ namespace Opm
             return;
 
         int phasePos;
+        Well2::GuideRateTarget wellTarget;
+        Group2::GuideRateTarget groupTarget;
+
+
         switch (injectorType) {
         case Well2::InjectorType::WATER:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Aqua];
+            wellTarget = Well2::GuideRateTarget::WAT;
+            groupTarget = Group2::GuideRateTarget::WAT;
             break;
         }
         case Well2::InjectorType::OIL:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Liquid];
+            wellTarget = Well2::GuideRateTarget::OIL;
+            groupTarget = Group2::GuideRateTarget::OIL;
             break;
         }
         case Well2::InjectorType::GAS:
         {
             phasePos = pu.phase_pos[BlackoilPhases::Vapour];
+            wellTarget = Well2::GuideRateTarget::GAS;
+            groupTarget = Group2::GuideRateTarget::GAS;
             break;
         }
         default:
             throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
         }
 
-        double groupTotalGuideRate = 0.0;
         double groupTargetReduction = 0.0;
-        wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTotalGuideRate, groupTargetReduction);
-        double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-        double fraction = guideRate / groupTotalGuideRate;
+        wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTargetReduction);
+        double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, wellTarget, /*isInjector*/true);
+        wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, groupTarget, /*isInjector*/true, fraction);
 
         switch(currentGroupControl) {
         case Group2::InjectionCMode::NONE:
@@ -1042,6 +1051,7 @@ namespace Opm
             const auto& parent = schedule.getGroup2( group.parent(), current_step_ );
             if (group.getTransferGroupEfficiencyFactor())
                 efficiencyFactor *= group.getGroupEfficiencyFactor();
+
             assembleGroupInjectionControl(parent, well_state, schedule, summaryState, injectorType, control_eq, efficiencyFactor, deferred_logger);
             break;
         }
@@ -1079,12 +1089,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::ORAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            int phasePos = pu.phase_pos[Oil];                    
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::OIL, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::OIL, /*isInjector*/false, fraction);
 
             const double rate_target = std::max(0.0, groupcontrols.oil_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
@@ -1094,12 +1103,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::WRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::WAT, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::WAT, /*isInjector*/false, fraction);
 
             const double rate_target = std::max(0.0, groupcontrols.water_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
@@ -1109,12 +1117,11 @@ namespace Opm
         }
         case Group2::ProductionCMode::GRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::GAS, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::GAS, /*isInjector*/false, fraction);
 
             const double rate_target = std::max(0.0, groupcontrols.gas_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::gasCompIdx));
@@ -1124,15 +1131,14 @@ namespace Opm
         }
         case Group2::ProductionCMode::LRAT:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
             phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            guideRate += wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::LIQ, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::LIQ, /*isInjector*/false, fraction);
+
             const double rate_target = std::max(0.0, groupcontrols.liquid_target / efficiencyFactor - groupTargetReduction);
             assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
 
@@ -1148,18 +1154,16 @@ namespace Opm
         }
         case Group2::ProductionCMode::RESV:
         {
-            double groupTotalGuideRate = 0.0;
             double groupTargetReduction = 0.0;
             int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            double guideRate = wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
             phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            guideRate += wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
             phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeTotalGuideRate(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTotalGuideRate, groupTargetReduction);
-            guideRate += wellGroupHelpers::computeGuideRate(well_state, well, phasePos);
-            double fraction = guideRate / groupTotalGuideRate;
+            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+
+            double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, current_step_, Base::guide_rate_, Well2::GuideRateTarget::RES, /*isInjector*/false);
+            wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, current_step_, Base::guide_rate_, Group2::GuideRateTarget::RES, /*isInjector*/false, fraction);
 
             EvalWell total_rate(numWellEq_ + numEq, 0.); // reservoir rate
             std::vector<double> convert_coeff(number_of_phases_, 1.0);

opm/simulators/wells/WellGroupHelpers.hpp

@@ -56,68 +56,6 @@ namespace Opm {
             accumulateGroupEfficiencyFactor(schedule.getGroup2(group.parent(), reportStepIdx), schedule, reportStepIdx, factor);
     }
 
-
-    inline double sumWellRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
-
-        double rate = 0.0;
-        for (const std::string& groupName : group.groups()) {
-            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
-            rate += groupTmp.getGroupEfficiencyFactor()*sumWellRates(groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
-        }
-        const auto& end = wellState.wellMap().end();
-        for (const std::string& wellName : group.wells()) {
-            const auto& it = wellState.wellMap().find( wellName );
-            if (it == end)  // the well is not found
-                continue;
-
-            int well_index = it->second[0];
-
-            const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
-            //only count producers or injectors
-            if ( (wellEcl.isProducer() && injector) ||  (wellEcl.isInjector() && !injector))
-                continue;
-
-            double factor = wellEcl.getEfficiencyFactor();
-            const auto wellrate_index = well_index * wellState.numPhases();
-            if (injector)
-                rate += factor * wellState.wellRates()[ wellrate_index + phasePos];
-            else
-                rate -= factor * wellState.wellRates()[ wellrate_index + phasePos];
-        }
-        return rate;
-    }
-
-    inline double sumWellResRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
-
-        double rate = 0.0;
-        for (const std::string& groupName : group.groups()) {
-            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
-            rate += groupTmp.getGroupEfficiencyFactor() * sumWellResRates(groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
-        }
-        const auto& end = wellState.wellMap().end();
-        for (const std::string& wellName : group.wells()) {
-            const auto& it = wellState.wellMap().find( wellName );
-            if (it == end)  // the well is not found
-                continue;
-
-            int well_index = it->second[0];
-
-            const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
-            //only count producers or injectors
-            if ( (wellEcl.isProducer() && injector) ||  (wellEcl.isInjector() && !injector))
-                continue;
-
-            double factor = wellEcl.getEfficiencyFactor();
-
-            const auto wellrate_index = well_index * wellState.numPhases();
-            if (injector)
-                rate += factor * wellState.wellReservoirRates()[ wellrate_index + phasePos];
-            else
-                rate -= factor * wellState.wellReservoirRates()[ wellrate_index + phasePos];
-        }
-        return rate;
-    }
-
     inline void setGroupControl(const Group2& group, const Schedule& schedule, const int reportStepIdx, const bool injector, WellStateFullyImplicitBlackoil& wellState) {
 
         for (const std::string& groupName : group.groups()) {
@@ -137,8 +75,13 @@ namespace Opm {
 
             int well_index = it->second[0];
             const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
+
+            if (wellEcl.getStatus() == Well2::Status::SHUT)
+                continue;
+
             if (!wellEcl.isAvailableForGroupControl())
                 continue;
+
             if (wellEcl.isProducer() && !injector)
                 wellState.currentProductionControls()[well_index] = Well2::ProducerCMode::GRUP;
 
@@ -147,36 +90,14 @@ namespace Opm {
         }
     }
 
-    inline double computeGuideRate(const WellStateFullyImplicitBlackoil& wellState, const Well2& well, const int phasePos )
-    {                
-        const auto& end = wellState.wellMap().end();
-        const auto& it = wellState.wellMap().find( well.name() );
-        if (it == end)  // the well is not found
-            return 0.0;
 
-        int well_index = it->second[0];
-
-         Opm::Well2::GuideRateTarget GuideRatePhase = well.getGuideRatePhase();
-#warning handle the case when phase != GuideRatePhase
-        if (GuideRatePhase == Opm::Well2::GuideRateTarget::UNDEFINED) {
-            const auto wellrate_index = well_index * wellState.numPhases();
-            double guideRate = wellState.wellPotentials()[wellrate_index + phasePos];
-            //if (guideRate == 0)
-            //    guideRate = 1; // if the well potential is not defined set it to 1
-#warning Add support for GUIDERAT
-            return guideRate*well.getGuideRateScalingFactor();
-        }
-        return well.getGuideRate()*well.getGuideRateScalingFactor();
-    }
-
-    inline void computeTotalGuideRate(const Group2& group, const WellStateFullyImplicitBlackoil& wellState, const Schedule& schedule, const int reportStepIdx, const int phasePos, const bool isInjector, double& groupTotalGuideRate, double& groupTargetReduction )
+    inline void computeGroupTargetReduction(const Group2& group, const WellStateFullyImplicitBlackoil& wellState, const Schedule& schedule, const int reportStepIdx, const int phasePos, const bool isInjector, double& groupTargetReduction )
     {
 
         for (const std::string& groupName : group.groups()) {
             const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
-            computeTotalGuideRate(groupTmp, wellState, schedule, reportStepIdx, phasePos, isInjector, groupTotalGuideRate, groupTargetReduction);
+            computeGroupTargetReduction(groupTmp, wellState, schedule, reportStepIdx, phasePos, isInjector, groupTargetReduction);
         }
-#warning return groupGuideRate if set, else sum it from below
         for (const std::string& wellName : group.wells()) {
             const auto& wellTmp = schedule.getWell2(wellName, reportStepIdx);
 
@@ -186,6 +107,9 @@ namespace Opm {
             if (wellTmp.isInjector() && !isInjector)
                  continue;
 
+            if (wellTmp.getStatus() == Well2::Status::SHUT)
+                continue;
+
             const auto& end = wellState.wellMap().end();
                 const auto& it = wellState.wellMap().find( wellName );
                 if (it == end)  // the well is not found
@@ -193,22 +117,137 @@ namespace Opm {
 
             int well_index = it->second[0];
             const auto wellrate_index = well_index * wellState.numPhases();
-            // only include wells under group control
+            // add contributino from wells not under group control
             if (isInjector) {
-                if (wellState.currentInjectionControls()[well_index] == Well2::InjectorCMode::GRUP)
-                    groupTotalGuideRate += computeGuideRate(wellState, wellTmp, phasePos);
-                else
+                if (wellState.currentInjectionControls()[well_index] != Well2::InjectorCMode::GRUP)
                     groupTargetReduction += wellState.wellRates()[wellrate_index + phasePos];
             } else {
-                if (wellState.currentProductionControls()[well_index] ==  Well2::ProducerCMode::GRUP)
-                    groupTotalGuideRate += computeGuideRate(wellState, wellTmp, phasePos);
-                else
+                if (wellState.currentProductionControls()[well_index] !=  Well2::ProducerCMode::GRUP)
                     groupTargetReduction -= wellState.wellRates()[wellrate_index + phasePos];
             }
         }
     }
 
 
+    inline double sumWellPhaseRates(const std::vector<double>& rates, const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos,
+                                    const bool injector) {
+
+        double rate = 0.0;
+        for (const std::string& groupName : group.groups()) {
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+            rate += groupTmp.getGroupEfficiencyFactor()*sumWellPhaseRates(rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
+        }
+        const auto& end = wellState.wellMap().end();
+        for (const std::string& wellName : group.wells()) {
+            const auto& it = wellState.wellMap().find( wellName );
+            if (it == end)  // the well is not found
+                continue;
+
+            int well_index = it->second[0];
+
+            const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
+            //only count producers or injectors
+            if ( (wellEcl.isProducer() && injector) ||  (wellEcl.isInjector() && !injector))
+                continue;
+
+            if (wellEcl.getStatus() == Well2::Status::SHUT)
+                continue;
+
+            double factor = wellEcl.getEfficiencyFactor();
+            const auto wellrate_index = well_index * wellState.numPhases();
+            if (injector)
+                rate += factor * rates[ wellrate_index + phasePos];
+            else
+                rate -= factor * rates[ wellrate_index + phasePos];
+        }
+        return rate;
+    }
+
+    inline double sumWellRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
+        return sumWellPhaseRates(wellState.wellRates(), group, schedule, wellState, reportStepIdx, phasePos, injector);
+    }
+
+    inline double sumWellResRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
+        return sumWellPhaseRates(wellState.wellReservoirRates(), group, schedule, wellState, reportStepIdx, phasePos, injector);
+    }
+
+
+    inline void updateGuideRateForGroups(const Group2& group, const Schedule& schedule, const PhaseUsage& pu, const int reportStepIdx, const double& simTime, GuideRate* guideRate, WellStateFullyImplicitBlackoil& wellState) {
+        for (const std::string& groupName : group.groups()) {
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+            updateGuideRateForGroups(groupTmp, schedule, pu, reportStepIdx, simTime, guideRate, wellState);
+        }
+        bool isInjector = group.isInjectionGroup();
+        bool isProducer = group.isProductionGroup();
+
+        if (!isInjector && !isProducer)
+            return;
+
+        double oilPot = 0.0;
+        if (pu.phase_used[BlackoilPhases::Liquid])
+            oilPot = sumWellPhaseRates(wellState.wellPotentials(), group, schedule, wellState, reportStepIdx, pu.phase_pos[BlackoilPhases::Liquid], isInjector);
+
+        double gasPot = 0.0;
+        if (pu.phase_used[BlackoilPhases::Vapour])
+            gasPot = sumWellPhaseRates(wellState.wellPotentials(), group, schedule, wellState, reportStepIdx, pu.phase_pos[BlackoilPhases::Vapour], isInjector);
+
+        double waterPot = 0.0;
+        if (pu.phase_used[BlackoilPhases::Aqua])
+            waterPot = sumWellPhaseRates(wellState.wellPotentials(), group, schedule, wellState, reportStepIdx, pu.phase_pos[BlackoilPhases::Aqua], isInjector);
+
+        guideRate->compute(group.name(), reportStepIdx, simTime, oilPot, gasPot, waterPot);
+    }
+
+    inline double wellFractionFromGuideRates(const Well2& well, const Schedule& schedule, const int reportStepIdx, const GuideRate* guideRate, const Well2::GuideRateTarget& wellTarget, const bool isInjector) {
+        double groupTotalGuideRate = 0.0;
+        const Group2& groupTmp = schedule.getGroup2(well.groupName(), reportStepIdx);
+        for (const std::string& wellName : groupTmp.wells()) {
+            const auto& wellTmp = schedule.getWell2(wellName, reportStepIdx);
+
+            if (wellTmp.isProducer() && isInjector)
+                 continue;
+
+            if (wellTmp.isInjector() && !isInjector)
+                 continue;
+
+            if (wellTmp.getStatus() == Well2::Status::SHUT)
+                continue;
+
+            groupTotalGuideRate += guideRate->get(wellName, wellTarget);
+        }
+        double wellGuideRate = guideRate->get(well.name(), wellTarget);
+        return wellGuideRate / groupTotalGuideRate;
+    }
+
+    inline double groupFractionFromGuideRates(const Group2& group, const Schedule& schedule, const int reportStepIdx, const GuideRate* guideRate, const Group2::GuideRateTarget& groupTarget, const bool isInjector) {
+        double groupTotalGuideRate = 0.0;
+        const Group2& groupParent = schedule.getGroup2(group.parent(), reportStepIdx);
+        for (const std::string& groupName : groupParent.groups()) {
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+
+            if ( (groupTmp.isProductionGroup() && !isInjector) ||
+                 (groupTmp.isInjectionGroup() && isInjector) )
+            {
+                groupTotalGuideRate += guideRate->get(groupName, groupTarget);
+            }
+        }
+        double groupGuideRate = guideRate->get(group.name(), groupTarget);
+        return groupGuideRate / groupTotalGuideRate;
+    }
+
+    inline void accumulateGroupFractions(const std::string& groupName, const std::string& controlGroupName, const Schedule& schedule, const int reportStepIdx, const GuideRate* guideRate, const Group2::GuideRateTarget& groupTarget, const bool isInjector, double fraction) {
+
+        const Group2& group = schedule.getGroup2(groupName, reportStepIdx);
+        if (groupName != controlGroupName) {
+            fraction *= groupFractionFromGuideRates(group, schedule, reportStepIdx, guideRate, groupTarget, isInjector);
+            accumulateGroupFractions(group.parent(), controlGroupName, schedule, reportStepIdx, guideRate, groupTarget, isInjector, fraction);
+        }
+
+        return;
+    }
+
+
+
     } // namespace wellGroupHelpers
 
 }

opm/simulators/wells/WellInterface.hpp

@@ -30,6 +30,7 @@
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well2.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
@@ -143,6 +144,8 @@ namespace Opm
 
         void setVFPProperties(const VFPProperties<VFPInjProperties,VFPProdProperties>* vfp_properties_arg);
 
+        void setGuideRate(const GuideRate* guide_rate_arg);
+
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
@@ -353,6 +356,8 @@ namespace Opm
 
         const VFPProperties<VFPInjProperties,VFPProdProperties>* vfp_properties_;
 
+        const GuideRate* guide_rate_;
+
         double gravity_;
 
         // For the conversion between the surface volume rate and resrevoir voidage rate

opm/simulators/wells/WellInterface_impl.hpp

@@ -176,8 +176,13 @@ namespace Opm
         vfp_properties_ = vfp_properties_arg;
     }
 
-
-
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    setGuideRate(const GuideRate* guide_rate_arg)
+    {
+        guide_rate_ = guide_rate_arg;
+    }
 
 
     template<typename TypeTag>