Merge pull request #2144 from totto82/nupcolFix2

Various group control fixes

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -677,6 +677,7 @@ namespace Opm {
                         const int np = numPhases();
                         for (int p = 0; p < np; ++p) {
                             well_state_.wellRates()[np * w + p] = 0.;
+                            well_state_.wellReservoirRates()[np * w + p] = 0.;
                         }
                         continue;
                     } else {
@@ -785,10 +786,25 @@ namespace Opm {
                 calculateExplicitQuantities(local_deferredLogger);
                 prepareTimeStep(local_deferredLogger);
             }
+            updateWellControls(local_deferredLogger, true);
+
             // only check group controls for iterationIdx smaller then nupcol
             const int reportStepIdx = ebosSimulator_.episodeIndex();
             const int nupcol = schedule().getNupcol(reportStepIdx);
-            updateWellControls(local_deferredLogger, iterationIdx < nupcol);
+            if (iterationIdx < nupcol) {
+                if( localWellsActive() ) {
+                    const Group2& fieldGroup = schedule().getGroup2("FIELD", reportStepIdx);
+                    std::vector<double> groupTargetReduction(numPhases(), 0.0);
+                    wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, well_state_, groupTargetReduction);
+                    std::vector<double> groupTargetReductionInj(numPhases(), 0.0);
+                    wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, well_state_, groupTargetReductionInj);
+                    std::vector<double> rein(numPhases(), 0.0);
+                    wellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, well_state_, rein);
+                    double resv = 0.0;
+                    wellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_, resv);
+                }
+            }
+
             // Set the well primary variables based on the value of well solutions
             initPrimaryVariablesEvaluation();
 
@@ -1663,7 +1679,8 @@ namespace Opm {
             if (group.has_control(Group2::InjectionCMode::REIN))
             {
                 double production_Rate = 0.0;
-                production_Rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
+                const Group2& groupRein = schedule().getGroup2(controls.reinj_group, reportStepIdx);
+                production_Rate += wellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
 
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
@@ -1673,16 +1690,18 @@ namespace Opm {
                 }                    }
             if (group.has_control(Group2::InjectionCMode::VREP))
             {
-                double voidage_Rate = 0.0;
-                voidage_Rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
-                voidage_Rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
-                voidage_Rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
+                double voidage_rate = 0.0;
+                const Group2& groupVoidage = schedule().getGroup2(controls.voidage_group, reportStepIdx);
+                voidage_rate += wellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
+                voidage_rate += wellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
+                voidage_rate += wellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
 
+                double total_rate = 0.0;
+                total_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true);
+                total_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true);
+                total_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true);
 
-                double current_rate = 0.0;
-                current_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
-
-                if (controls.target_void_fraction*voidage_Rate < current_rate) {
+                if (controls.target_void_fraction*voidage_rate < total_rate) {
                     actionOnBrokenConstraints(group, Group2::InjectionCMode::VREP, reportStepIdx, deferred_logger);
                 }
             }

opm/simulators/wells/MultisegmentWell.hpp

@@ -418,6 +418,17 @@ namespace Opm
 
         void updateUpwindingSegments();
 
+        // turn on crossflow to avoid singular well equations
+        // when the well is banned from cross-flow and the BHP is not properly initialized,
+        // we turn on crossflow to avoid singular well equations. It can result in wrong-signed
+        // well rates, it can cause problem for THP calculation
+        // TODO: looking for better alternative to avoid wrong-signed well rates
+        bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
+
+        // for a well, when all drawdown are in the wrong direction, then this well will not
+        // be able to produce/inject .
+        bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
+
     };
 
 }

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1980,7 +1980,17 @@ namespace Opm
         const auto& well = well_ecl_;
         const auto pu = phaseUsage();
         const Group2::InjectionCMode& currentGroupControl = well_state.currentInjectionGroupControl(group.name());
-        if (!group.isInjectionGroup() && currentGroupControl!=Group2::InjectionCMode::FLD )
+        if (currentGroupControl == Group2::InjectionCMode::FLD) {
+            // Inject share of parents control
+            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);
+            return;
+        }
+
+        if (!group.isInjectionGroup())
             return;
 
         const auto& groupcontrols = group.injectionControls(summaryState);
@@ -2019,8 +2029,8 @@ namespace Opm
             throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
         }
 
-        double groupTargetReduction = 0.0;
-        wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTargetReduction);
+        const std::vector<double>& groupTargetReductions = well_state.currentProductionGroupReductionRates(group.name());
+        double groupTargetReduction = groupTargetReductions[phasePos];
         double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, wellTarget, /*isInjector*/true);
         wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, groupTarget, /*isInjector*/true, fraction);
 
@@ -2047,9 +2057,7 @@ namespace Opm
         }
         case Group2::InjectionCMode::REIN:
         {
-            double productionRate = 0.0;
-            const Group2& reinGroup = schedule.getGroup2( groupcontrols.reinj_group, current_step_ );
-            productionRate += wellGroupHelpers::sumWellRates(reinGroup, schedule, well_state, current_step_, phasePos, /*isInjector*/false);
+            double productionRate = well_state.currentInjectionVREPRates(groupcontrols.reinj_group);
             productionRate /= efficiencyFactor;
             double target = std::max(0.0, (groupcontrols.target_reinj_fraction*productionRate - groupTargetReduction));
             control_eq = getSegmentGTotal(0) / scaling - fraction * target;
@@ -2060,11 +2068,7 @@ namespace Opm
             std::vector<double> convert_coeff(number_of_phases_, 1.0);
             Base::rateConverter_.calcCoeff(/*fipreg*/ 0, Base::pvtRegionIdx_, convert_coeff);
             double coeff = convert_coeff[phasePos];
-            double voidageRate = 0.0;
-            const Group2& voidageGroup = schedule.getGroup2( groupcontrols.voidage_group, current_step_ );
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Aqua], /*injector*/false);
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Liquid], /*injector*/false);
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Vapour], /*injector*/false);
+            double voidageRate = well_state.currentInjectionVREPRates(groupcontrols.voidage_group);
             voidageRate /= efficiencyFactor;
             double target = std::max(0.0, ( groupcontrols.target_void_fraction*voidageRate/coeff - groupTargetReduction));
             control_eq = getSegmentGTotal(0) / scaling - fraction * target ;
@@ -2072,11 +2076,8 @@ namespace Opm
         }
         case Group2::InjectionCMode::FLD:
         {
-            // Inject share of parents control
-            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);
+            // The FLD case is handled earlier
+            assert(false);
             break;
         }
 
@@ -2100,10 +2101,21 @@ namespace Opm
         const auto pu = phaseUsage();
 
         const Group2::ProductionCMode& currentGroupControl = well_state.currentProductionGroupControl(group.name());
-        if (!group.isProductionGroup() && currentGroupControl != Group2::ProductionCMode::FLD)
+
+        if (currentGroupControl == Group2::ProductionCMode::FLD ) {
+            // Produce share of parents control
+            const auto& parent = schedule.getGroup2( group.parent(), current_step_ );
+            if (group.getTransferGroupEfficiencyFactor())
+                efficiencyFactor *= group.getGroupEfficiencyFactor();
+
+            assembleGroupProductionControl(parent, well_state, schedule, summaryState, control_eq, efficiencyFactor, deferred_logger);
+            return;
+        }
+        if (!group.isProductionGroup())
             return;
 
         const auto& groupcontrols = group.productionControls(summaryState);
+        const std::vector<double>& groupTargetReductions = well_state.currentProductionGroupReductionRates(group.name());
 
         switch(currentGroupControl) {
         case Group2::ProductionCMode::NONE:
@@ -2113,9 +2125,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::ORAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Oil]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::OIL, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::OIL, /*isInjector*/false, fraction);
 
@@ -2127,9 +2137,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::WRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Water]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::WAT, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::WAT, /*isInjector*/false, fraction);
 
@@ -2141,9 +2149,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::GRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Gas]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::GAS, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::GAS, /*isInjector*/false, fraction);
             const double rate_target = std::max(0.0, groupcontrols.gas_target / efficiencyFactor - groupTargetReduction);
@@ -2154,11 +2160,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::LRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
-            phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Oil]] + groupTargetReductions[pu.phase_pos[Water]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::LIQ, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::LIQ, /*isInjector*/false, fraction);
 
@@ -2187,15 +2189,10 @@ namespace Opm
         }
         case Group2::ProductionCMode::FLD:
         {
-            // Produce share of parents control
-            const auto& parent = schedule.getGroup2( group.parent(), current_step_ );
-            if (group.getTransferGroupEfficiencyFactor())
-                efficiencyFactor *= group.getGroupEfficiencyFactor();
-
-            assembleGroupProductionControl(parent, well_state, schedule, summaryState, control_eq, efficiencyFactor, deferred_logger);
+            // The FLD case is handled earlier
+            assert(false);
             break;
         }
-
         default:
             OPM_DEFLOG_THROW(std::runtime_error, "Unvallied group control specified for group "  + well.groupName(), deferred_logger );
         }
@@ -2609,7 +2606,7 @@ namespace Opm
         //
         // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
 
-        const bool allow_cf = getAllowCrossFlow();
+        const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
 
         const int nseg = numberOfSegments();
 
@@ -2726,9 +2723,55 @@ namespace Opm
             }
         }
     }
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
+    {
+        return !getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
+    }
 
 
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    allDrawDownWrongDirection(const Simulator& ebos_simulator) const
+    {
+        bool all_drawdown_wrong_direction = true;
+        const int nseg = numberOfSegments();
 
+        for (int seg = 0; seg < nseg; ++seg) {
+            const EvalWell segment_pressure = getSegmentPressure(seg);
+            for (const int perf : segment_perforations_[seg]) {
+
+                const int cell_idx = well_cells_[perf];
+                const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
+                const auto& fs = intQuants.fluidState();
+
+                // pressure difference between the segment and the perforation
+                const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
+                // pressure difference between the perforation and the grid cell
+                const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];
+
+                const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
+                const double perf_press = pressure_cell - cell_perf_press_diff;
+                // Pressure drawdown (also used to determine direction of flow)
+                // TODO: not 100% sure about the sign of the seg_perf_press_diff
+                const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
+
+                // for now, if there is one perforation can produce/inject in the correct
+                // direction, we consider this well can still produce/inject.
+                // TODO: it can be more complicated than this to cause wrong-signed rates
+                if ( (drawdown < 0. && well_type_ == INJECTOR) ||
+                     (drawdown > 0. && well_type_ == PRODUCER) )  {
+                    all_drawdown_wrong_direction = false;
+                    break;
+                }
+            }
+        }
+
+        return all_drawdown_wrong_direction;
+    }
 
 
     template<typename TypeTag>

opm/simulators/wells/StandardWell_impl.hpp

@@ -961,15 +961,24 @@ namespace Opm
         const auto pu = phaseUsage();
         const auto& groupcontrols = group.injectionControls(summaryState);
         const Group2::InjectionCMode& currentGroupControl = well_state.currentInjectionGroupControl(group.name());
-        // it may have a injection group higher up in tree
-        if (!group.isInjectionGroup() && currentGroupControl != Group2::InjectionCMode::FLD)
+
+        if (currentGroupControl == Group2::InjectionCMode::FLD) {
+            // Inject share of parents control
+            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);
+            return;
+        }
+
+        if (!group.isInjectionGroup())
             return;
 
         int phasePos;
         Well2::GuideRateTarget wellTarget;
         Group2::GuideRateTarget groupTarget;
 
-
         switch (injectorType) {
         case Well2::InjectorType::WATER:
         {
@@ -996,8 +1005,8 @@ namespace Opm
             throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
         }
 
-        double groupTargetReduction = 0.0;
-        wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/true, groupTargetReduction);
+        const std::vector<double>& groupInjectionReductions = well_state.currentInjectionGroupReductionRates(group.name());
+        double groupTargetReduction = groupInjectionReductions[phasePos];
         double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, wellTarget, /*isInjector*/true);
         wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, groupTarget, /*isInjector*/true, fraction);
 
@@ -1024,9 +1033,7 @@ namespace Opm
         }
         case Group2::InjectionCMode::REIN:
         {
-            double productionRate = 0.0;
-            const Group2& reinGroup = schedule.getGroup2( groupcontrols.reinj_group, current_step_ );
-            productionRate += wellGroupHelpers::sumWellRates(reinGroup, schedule, well_state, current_step_, phasePos, /*isInjector*/false);
+            double productionRate = well_state.currentInjectionVREPRates(groupcontrols.reinj_group);
             productionRate /= efficiencyFactor;
             double target = std::max(0.0, (groupcontrols.target_reinj_fraction*productionRate - groupTargetReduction));
             control_eq = getWQTotal() - fraction * target;            
@@ -1037,27 +1044,33 @@ namespace Opm
             std::vector<double> convert_coeff(number_of_phases_, 1.0);
             Base::rateConverter_.calcCoeff(/*fipreg*/ 0, Base::pvtRegionIdx_, convert_coeff);
             double coeff = convert_coeff[phasePos];
-            double voidageRate = 0.0;
-            const Group2& voidageGroup = schedule.getGroup2( groupcontrols.voidage_group, current_step_ );
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Aqua], /*injector*/false);
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Liquid], /*injector*/false);
-            voidageRate += wellGroupHelpers::sumWellResRates(voidageGroup, schedule, well_state, current_step_, pu.phase_pos[BlackoilPhases::Vapour], /*injector*/false);
+            double voidageRate = well_state.currentInjectionVREPRates(groupcontrols.voidage_group);
+
+            double injReduction = 0.0;
+
+            if (groupcontrols.phase != Phase::WATER)
+                injReduction += groupInjectionReductions[pu.phase_pos[BlackoilPhases::Aqua]]*convert_coeff[pu.phase_pos[BlackoilPhases::Aqua]];
+
+            if (groupcontrols.phase != Phase::OIL)
+                injReduction += groupInjectionReductions[pu.phase_pos[BlackoilPhases::Liquid]]*convert_coeff[pu.phase_pos[BlackoilPhases::Liquid]];
+
+            if (groupcontrols.phase != Phase::GAS)
+                injReduction += groupInjectionReductions[pu.phase_pos[BlackoilPhases::Vapour]]*convert_coeff[pu.phase_pos[BlackoilPhases::Vapour]];
+
+            voidageRate -= injReduction;
+
             voidageRate /= efficiencyFactor;
+
             double target = std::max(0.0, ( groupcontrols.target_void_fraction*voidageRate/coeff - groupTargetReduction));
             control_eq = getWQTotal() - fraction * target;
             break;
         }
         case Group2::InjectionCMode::FLD:
         {
-            // Inject share of parents control
-            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);
+            // The FLD case is handled earlier
+            assert(false);
             break;
         }
-
         default:
             OPM_DEFLOG_THROW(std::runtime_error, "Unvalid group control specified for group "  + well.groupName(), deferred_logger );
         }
@@ -1078,11 +1091,23 @@ namespace Opm
         const auto pu = phaseUsage();
 
         const Group2::ProductionCMode& currentGroupControl = well_state.currentProductionGroupControl(group.name());
-        if (!group.isProductionGroup() && currentGroupControl != Group2::ProductionCMode::FLD)
+        if (currentGroupControl == Group2::ProductionCMode::FLD) {
+            // Produce share of parents control
+            const auto& parent = schedule.getGroup2( group.parent(), current_step_ );
+            if (group.getTransferGroupEfficiencyFactor())
+                efficiencyFactor *= group.getGroupEfficiencyFactor();
+
+            assembleGroupProductionControl(parent, well_state, schedule, summaryState, control_eq, efficiencyFactor, deferred_logger);
+            return;
+        }
+
+        if (!group.isProductionGroup())
             return;
 
         const auto& groupcontrols = group.productionControls(summaryState);
 
+        const std::vector<double>& groupTargetReductions = well_state.currentProductionGroupReductionRates(group.name());
+
         switch(currentGroupControl) {
         case Group2::ProductionCMode::NONE:
         {
@@ -1091,9 +1116,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::ORAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Oil];                    
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Oil]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::OIL, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::OIL, /*isInjector*/false, fraction);
 
@@ -1105,9 +1128,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::WRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Water]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::WAT, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::WAT, /*isInjector*/false, fraction);
 
@@ -1119,9 +1140,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::GRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Gas]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::GAS, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::GAS, /*isInjector*/false, fraction);
 
@@ -1133,11 +1152,7 @@ namespace Opm
         }
         case Group2::ProductionCMode::LRAT:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
-            phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Oil]] + groupTargetReductions[pu.phase_pos[Water]];
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::LIQ, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::LIQ, /*isInjector*/false, fraction);
 
@@ -1156,13 +1171,9 @@ namespace Opm
         }
         case Group2::ProductionCMode::RESV:
         {
-            double groupTargetReduction = 0.0;
-            int phasePos = pu.phase_pos[Water];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
-            phasePos = pu.phase_pos[Oil];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
-            phasePos = pu.phase_pos[Gas];
-            wellGroupHelpers::computeGroupTargetReduction(group, well_state, schedule, current_step_, phasePos, /*isInjector*/false, groupTargetReduction);
+            double groupTargetReduction = groupTargetReductions[pu.phase_pos[Oil]]
+                    + groupTargetReductions[pu.phase_pos[Gas]]
+                    + groupTargetReductions[pu.phase_pos[Water]];
 
             double fraction = wellGroupHelpers::wellFractionFromGuideRates(well, schedule, well_state, current_step_, Base::guide_rate_, Well2::GuideRateTarget::RES, /*isInjector*/false);
             wellGroupHelpers::accumulateGroupFractions(well.groupName(), group.name(), schedule, well_state, current_step_, Base::guide_rate_, Group2::GuideRateTarget::RES, /*isInjector*/false, fraction);
@@ -1186,12 +1197,8 @@ namespace Opm
         }
         case Group2::ProductionCMode::FLD:
         {
-            // Produce share of parents control
-            const auto& parent = schedule.getGroup2( group.parent(), current_step_ );
-            if (group.getTransferGroupEfficiencyFactor())
-                efficiencyFactor *= group.getGroupEfficiencyFactor();
-
-            assembleGroupProductionControl(parent, well_state, schedule, summaryState, control_eq, efficiencyFactor, deferred_logger);
+            // The FLD case is handled earlier
+            assert(false);
             break;
         }
 

opm/simulators/wells/WellGroupHelpers.hpp

@@ -35,17 +35,21 @@ namespace Opm {
             setCmodeGroup( schedule.getGroup2(groupName, reportStepIdx), schedule, summaryState, reportStepIdx, wellState);
         }
 
-        // use FLD as default control
-        wellState.setCurrentInjectionGroupControl(group.name(), Group2::InjectionCMode::FLD);
-        wellState.setCurrentProductionGroupControl(group.name(), Group2::ProductionCMode::FLD);
-
-        if (group.isInjectionGroup()) {
-            const auto controls = group.injectionControls(summaryState);
-            wellState.setCurrentInjectionGroupControl(group.name(), controls.cmode);
+        if (!wellState.hasInjectionGroupControl(group.name())) {
+            // use NONE as default control
+            wellState.setCurrentInjectionGroupControl(group.name(), Group2::InjectionCMode::NONE);
+            if (group.isInjectionGroup()) {
+                const auto controls = group.injectionControls(summaryState);
+                wellState.setCurrentInjectionGroupControl(group.name(), controls.cmode);
+            }
         }
-        if (group.isProductionGroup()) {
-            const auto controls = group.productionControls(summaryState);
-            wellState.setCurrentProductionGroupControl(group.name(), controls.cmode);
+        if (!wellState.hasProductionGroupControl(group.name())) {
+            // use NONE as default control
+            wellState.setCurrentProductionGroupControl(group.name(), Group2::ProductionCMode::NONE);
+            if (group.isProductionGroup()) {
+                const auto controls = group.productionControls(summaryState);
+                wellState.setCurrentProductionGroupControl(group.name(), controls.cmode);
+            }
         }
     }
 
@@ -90,6 +94,10 @@ namespace Opm {
             }
 
             if (wellEcl.isInjector() && injector) {
+                // only switch if the well phase is the same as the group phase
+                if (group.injection_phase() != wellEcl.getPreferredPhase())
+                    continue;
+
                 if (wellState.currentInjectionControls()[well_index] != Well2::InjectorCMode::GRUP) {
                     wellState.currentInjectionControls()[well_index] = Well2::InjectorCMode::GRUP;
                     ss <<"\n Injector " << wellName << " switches to GRUP control limit";
@@ -211,6 +219,76 @@ namespace Opm {
         return rate;
     }
 
+    inline void updateGroupTargetReduction(const Group2& group, const Schedule& schedule, const int reportStepIdx, const bool isInjector, WellStateFullyImplicitBlackoil& wellState, std::vector<double>& groupTargetReduction)
+    {
+        const int np = wellState.numPhases();
+        for (const std::string& groupName : group.groups()) {
+            std::vector<double> thisGroupTargetReduction(np, 0.0);
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+            updateGroupTargetReduction(groupTmp, schedule, reportStepIdx, isInjector, wellState, thisGroupTargetReduction);
+
+            // accumulate group contribution from sub group
+            if (isInjector) {
+                const Group2::InjectionCMode& currentGroupControl = wellState.currentInjectionGroupControl(groupName);
+                if (currentGroupControl != Group2::InjectionCMode::FLD) {
+                    for (int phase = 0; phase < np; phase++) {
+                        groupTargetReduction[phase] += sumWellRates(groupTmp, schedule, wellState, reportStepIdx, phase, isInjector);
+                    }
+                    continue;
+                }
+            } else {
+                const Group2::ProductionCMode& currentGroupControl = wellState.currentProductionGroupControl(groupName);
+                if (currentGroupControl != Group2::ProductionCMode::FLD) {
+                    for (int phase = 0; phase < np; phase++) {
+                        groupTargetReduction[phase] += sumWellRates(groupTmp, schedule, wellState, reportStepIdx, phase, isInjector);
+                    }
+                    continue;
+                }
+            }
+            // or accumulate directly from the wells if controled from its parents
+            for (int phase = 0; phase < np; phase++) {
+                groupTargetReduction[phase] += thisGroupTargetReduction[phase];
+            }
+        }
+        for (const std::string& wellName : group.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;
+
+            const auto& end = wellState.wellMap().end();
+            const auto& it = wellState.wellMap().find( wellName );
+            if (it == end)  // the well is not found
+                continue;
+
+            int well_index = it->second[0];
+            const auto wellrate_index = well_index * wellState.numPhases();
+            // add contributino from wells not under group control
+            if (isInjector) {
+                if (wellState.currentInjectionControls()[well_index] != Well2::InjectorCMode::GRUP)
+                    for (int phase = 0; phase < np; phase++) {
+                        groupTargetReduction[phase] += wellState.wellRates()[wellrate_index + phase];
+                    }
+            } else {
+                if (wellState.currentProductionControls()[well_index] !=  Well2::ProducerCMode::GRUP)
+                    for (int phase = 0; phase < np; phase++) {
+                        groupTargetReduction[phase] -= wellState.wellRates()[wellrate_index + phase];
+                    }
+            }
+        }
+        if (isInjector)
+            wellState.setCurrentInjectionGroupReductionRates(group.name(), groupTargetReduction);
+        else
+            wellState.setCurrentProductionGroupReductionRates(group.name(), groupTargetReduction);
+    }
+
+
 
     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()) {
@@ -238,6 +316,37 @@ namespace Opm {
         guideRate->compute(group.name(), reportStepIdx, simTime, oilPot, gasPot, waterPot);
     }
 
+    inline void updateVREPForGroups(const Group2& group, const Schedule& schedule, const int reportStepIdx, WellStateFullyImplicitBlackoil& wellState, double& resv) {
+        for (const std::string& groupName : group.groups()) {
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+            double thisResv = 0.0;
+            updateVREPForGroups(groupTmp, schedule, reportStepIdx, wellState, thisResv);
+            resv += thisResv;
+        }
+        const int np = wellState.numPhases();
+        for (int phase = 0; phase < np; ++phase) {
+            resv += sumWellPhaseRates(wellState.wellReservoirRates(), group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false);
+        }
+
+        wellState.setCurrentInjectionVREPRates(group.name(), resv);
+    }
+
+    inline void updateREINForGroups(const Group2& group, const Schedule& schedule, const int reportStepIdx, WellStateFullyImplicitBlackoil& wellState, std::vector<double>& rein) {
+        const int np = wellState.numPhases();
+        for (const std::string& groupName : group.groups()) {
+            const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
+            std::vector<double> thisRein(np, 0.0);
+            updateREINForGroups(groupTmp, schedule, reportStepIdx, wellState, thisRein);
+            for (int phase = 0; phase < np; ++phase) {
+                rein[phase] = thisRein[phase];
+            }
+        }
+        for (int phase = 0; phase < np; ++phase) {
+            rein[phase] = sumWellPhaseRates(wellState.wellRates(), group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false);
+        }
+        wellState.setCurrentInjectionREINRates(group.name(), rein);
+    }
+
     inline double wellFractionFromGuideRates(const Well2& well, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, 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);

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -160,9 +160,6 @@ namespace Opm
 
             current_injection_controls_.resize(nw);
             current_production_controls_.resize(nw);
-            current_injection_group_controls_.clear();
-            current_production_group_controls_.clear();
-
 
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
@@ -206,6 +203,12 @@ namespace Opm
                             wellRates()[ idx ] = prevState->wellRates()[ oldidx ];
                         }
 
+                        // wellResrates
+                        for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
+                        {
+                            wellReservoirRates()[ idx ] = prevState->wellReservoirRates()[ oldidx ];
+                        }
+
                         // perfPhaseRates
                         const int oldPerf_idx_beg = (*it).second[ 1 ];
                         const int num_perf_old_well = (*it).second[ 2 ];
@@ -326,6 +329,14 @@ namespace Opm
         std::vector<Well2::ProducerCMode>& currentProductionControls() { return current_production_controls_; }
         const std::vector<Well2::ProducerCMode>& currentProductionControls() const { return current_production_controls_; }
 
+        bool hasProductionGroupControl(const std::string& groupName) {
+            return current_production_group_controls_.count(groupName) > 0;
+        }
+
+        bool hasInjectionGroupControl(const std::string& groupName) {
+            return current_injection_group_controls_.count(groupName) > 0;
+        }
+
         /// One current control per group.
         void setCurrentProductionGroupControl(const std::string& groupName, const Group2::ProductionCMode& groupControl ) {
             current_production_group_controls_[groupName] = groupControl;
@@ -334,7 +345,7 @@ namespace Opm
             auto it = current_production_group_controls_.find(groupName);
 
             if (it == current_production_group_controls_.end())
-                OPM_THROW(std::logic_error, "Could not find any well control for " << groupName);
+                OPM_THROW(std::logic_error, "Could not find any control for production group " << groupName);
 
             return it->second;
         }
@@ -347,11 +358,61 @@ namespace Opm
             auto it = current_injection_group_controls_.find(groupName);
 
             if (it == current_injection_group_controls_.end())
-                OPM_THROW(std::logic_error, "Could not find any well control for " << groupName);
+                OPM_THROW(std::logic_error, "Could not find any control for injection group " << groupName);
 
             return it->second;
         }
 
+
+        void setCurrentProductionGroupReductionRates(const std::string& groupName, const std::vector<double>& target ) {
+            production_group_reduction_rates[groupName] = target;
+        }
+        const std::vector<double>& currentProductionGroupReductionRates(const std::string& groupName) const {
+            auto it = production_group_reduction_rates.find(groupName);
+
+            if (it == production_group_reduction_rates.end())
+                OPM_THROW(std::logic_error, "Could not find any reduction rates for production group  " << groupName);
+
+            return it->second;
+        }
+
+        void setCurrentInjectionGroupReductionRates(const std::string& groupName, const std::vector<double>& target ) {
+            injection_group_reduction_rates[groupName] = target;
+        }
+        const std::vector<double>& currentInjectionGroupReductionRates(const std::string& groupName) const {
+            auto it = injection_group_reduction_rates.find(groupName);
+
+            if (it == injection_group_reduction_rates.end())
+                OPM_THROW(std::logic_error, "Could not find any reduction rates for injection group " << groupName);
+
+            return it->second;
+        }
+
+        void setCurrentInjectionVREPRates(const std::string& groupName, const double& target ) {
+            injection_group_vrep_rates[groupName] = target;
+        }
+        const double& currentInjectionVREPRates(const std::string& groupName) const {
+            auto it = injection_group_vrep_rates.find(groupName);
+
+            if (it == injection_group_vrep_rates.end())
+                OPM_THROW(std::logic_error, "Could not find any VREP rates for group " << groupName);
+
+            return it->second;
+        }
+
+        void setCurrentInjectionREINRates(const std::string& groupName, const std::vector<double>& target ) {
+            injection_group_rein_rates[groupName] = target;
+        }
+        const std::vector<double>& currentInjectionREINRates(const std::string& groupName) const {
+            auto it = injection_group_rein_rates.find(groupName);
+
+            if (it == injection_group_rein_rates.end())
+                OPM_THROW(std::logic_error, "Could not find any REIN rates for group " << groupName);
+
+            return it->second;
+        }
+
+
         data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override
         {
             data::Wells res = WellState::report(pu, globalCellIdxMap);
@@ -783,6 +844,10 @@ namespace Opm
         std::map<std::string, Group2::ProductionCMode> current_production_group_controls_;
         std::map<std::string, Group2::InjectionCMode> current_injection_group_controls_;
 
+        std::map<std::string, std::vector<double>> production_group_reduction_rates;
+        std::map<std::string, std::vector<double>> injection_group_reduction_rates;
+        std::map<std::string, double> injection_group_vrep_rates;
+        std::map<std::string, std::vector<double>> injection_group_rein_rates;
 
         std::vector<double> perfRateSolvent_;