Merge pull request #2223 from totto82/makeGrupMPI

communicate group rates

compareECLFiles.cmake

@@ -469,6 +469,12 @@ if(MPI_FOUND)
                                        ABS_TOL ${abs_tol_parallel}
                                        REL_TOL ${rel_tol_parallel})
 
+  add_test_compare_parallel_simulation(CASENAME spe9group
+                                       FILENAME SPE9_CP_GROUP
+                                       SIMULATOR flow
+                                       ABS_TOL ${abs_tol_parallel}
+                                       REL_TOL ${coarse_rel_tol_parallel})
+
   add_test_compare_parallel_simulation(CASENAME spe3
                                        FILENAME SPE3CASE1
                                        SIMULATOR flow

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -209,6 +209,7 @@ namespace Opm {
 
         const Grid& grid = ebosSimulator_.vanguard().grid();
         const auto& summaryState = ebosSimulator_.vanguard().summaryState();
+        int globalNumWells = 0;
         // Make wells_ecl_ contain only this partition's non-shut wells.
         {
             const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
@@ -216,6 +217,7 @@ namespace Opm {
                 return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
             };
             auto w = schedule().getWells(timeStepIdx);
+            globalNumWells = w.size();
             w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
             wells_ecl_.swap(w);
         }
@@ -261,7 +263,7 @@ namespace Opm {
             }
             cellPressures[cellIdx] = perf_pressure;
         }
-        well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState);
+        well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState, globalNumWells);
 
         // handling MS well related
         if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
@@ -377,20 +379,8 @@ namespace Opm {
         }
 
         //compute well guideRates
-        const int np = numPhases();
-        for (const auto& well : well_container_) {
-            const auto wpot = well_state_.wellPotentials().data() + well->indexOfWell()*np;
-            const double oilpot = (phase_usage_.phase_used[BlackoilPhases::Liquid] > 0)
-                ? wpot[phase_usage_.phase_pos[BlackoilPhases::Liquid]]
-                : 0.0;
-            const double gaspot = (phase_usage_.phase_used[BlackoilPhases::Vapour] > 0)
-                ? wpot[phase_usage_.phase_pos[BlackoilPhases::Vapour]]
-                : 0.0;
-            const double waterpot = (phase_usage_.phase_used[BlackoilPhases::Aqua] > 0)
-                ? wpot[phase_usage_.phase_pos[BlackoilPhases::Aqua]]
-                : 0.0;
-            guideRate_->compute(well->name(), reportStepIdx, simulationTime, oilpot, gaspot, waterpot);
-        }
+        const auto& comm = ebosSimulator_.vanguard().grid().comm();
+        wellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, well_state_, comm, guideRate_.get());      
     }
 
 
@@ -525,7 +515,7 @@ namespace Opm {
         // wells that have information written to the restart file.
         const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0);
         const auto& summaryState = ebosSimulator_.vanguard().summaryState();
-
+        int globalNumWells = 0;
         // Make wells_ecl_ contain only this partition's non-shut wells.
         {
             const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
@@ -533,6 +523,7 @@ namespace Opm {
                 return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
             };
             auto w = schedule().getWells(report_step);
+            globalNumWells = w.size();
             w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
             wells_ecl_.swap(w);
         }
@@ -544,13 +535,13 @@ namespace Opm {
             const auto phaseUsage = phaseUsageFromDeck(eclState());
             const size_t numCells = Opm::UgGridHelpers::numCells(grid());
             const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal());
-            well_state_.resize(wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState); // Resize for restart step
+            well_state_.resize(wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState, globalNumWells); // Resize for restart step
             wellsToState(restartValues.wells, phaseUsage, handle_ms_well, well_state_);
         }
 
         // for ecl compatible restart the current controls are not written
         const auto& ioCfg = eclState().getIOConfig();
-        const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST();
+        const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST();        
         if (true || ecl_compatible_rst) { // always set the control from the schedule
             for (int w = 0; w <nw; ++w) {
                 const auto& well = wells_ecl_[w];
@@ -1185,26 +1176,27 @@ namespace Opm {
 
         // the group target reduction rates needs to be update since wells may have swicthed to/from GRUP control
         // Currently the group targer reduction does not honor NUPCOL
-        if( localWellsActive() ) {
-            std::vector<double> groupTargetReduction(numPhases(), 0.0);
-            wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, well_state_nupcol_, well_state_, groupTargetReduction);
-            std::vector<double> groupTargetReductionInj(numPhases(), 0.0);
-            wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, well_state_nupcol_, well_state_, groupTargetReductionInj);
-        }
+        std::vector<double> groupTargetReduction(numPhases(), 0.0);
+        wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, well_state_nupcol_, well_state_, groupTargetReduction);
+        std::vector<double> groupTargetReductionInj(numPhases(), 0.0);
+        wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, well_state_nupcol_, well_state_, groupTargetReductionInj);
+
+        const auto& comm = ebosSimulator_.vanguard().grid().comm();
+        well_state_.communicateGroupReductionRates(comm);
+        well_state_.updateGlobalIsGrup(schedule(), reportStepIdx, comm);
 
         const double simulationTime = ebosSimulator_.time();
         std::vector<double> pot(numPhases(), 0.0);
-        wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, /*isInjector*/ false, well_state_, guideRate_.get(), pot);
+        wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, /*isInjector*/ false, well_state_, comm, guideRate_.get(), pot);
         std::vector<double> potInj(numPhases(), 0.0);
-        wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, /*isInjector*/ true, well_state_, guideRate_.get(), potInj);
+        wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, /*isInjector*/ true, well_state_, comm, guideRate_.get(), potInj);
 
-        if( localWellsActive() ) {
-            std::vector<double> rein(numPhases(), 0.0);
-            const auto& summaryState = ebosSimulator_.vanguard().summaryState();
-            wellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, phase_usage_, summaryState, well_state_nupcol_, well_state_, rein);
-            double resv = 0.0;
-            wellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol_, well_state_, resv);
-        }
+        std::vector<double> rein(numPhases(), 0.0);
+        const auto& summaryState = ebosSimulator_.vanguard().summaryState();
+        wellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, phase_usage_, summaryState, well_state_nupcol_, well_state_, rein);
+        double resv = 0.0;
+        wellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol_, well_state_, resv);
+        well_state_.communicateReinVrep(comm);
 
         // compute wsolvent fraction for REIN wells
         updateWsolvent(fieldGroup, schedule(), reportStepIdx,  well_state_nupcol_);
@@ -1690,6 +1682,7 @@ namespace Opm {
 
         const int reportStepIdx = ebosSimulator_.episodeIndex();
         const auto& summaryState = ebosSimulator_.vanguard().summaryState();
+        const auto& comm = ebosSimulator_.vanguard().grid().comm();
         auto& well_state = well_state_;
 
         if (group.isInjectionGroup())
@@ -1725,6 +1718,10 @@ namespace Opm {
             {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
+
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.surface_max_rate < current_rate) {
                     actionOnBrokenConstraints(group, Group::InjectionCMode::RATE, reportStepIdx, deferred_logger);
                 }
@@ -1733,6 +1730,9 @@ namespace Opm {
             {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.resv_max_rate < current_rate) {
                     actionOnBrokenConstraints(group, Group::InjectionCMode::RESV, reportStepIdx, deferred_logger);
                 }                    }
@@ -1742,9 +1742,15 @@ namespace Opm {
                 const Group& groupRein = schedule().getGroup(controls.reinj_group, reportStepIdx);
                 production_Rate += wellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
 
+                // sum over all nodes
+                production_Rate = comm.sum(production_Rate);
+
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
 
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.target_reinj_fraction*production_Rate < current_rate) {
                     actionOnBrokenConstraints(group, Group::InjectionCMode::REIN, reportStepIdx, deferred_logger);
                 }                    }
@@ -1756,11 +1762,17 @@ namespace Opm {
                 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);
 
+                // sum over all nodes
+                voidage_rate = comm.sum(voidage_rate);
+
                 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);
 
+                // sum over all nodes
+                total_rate = comm.sum(total_rate);
+
                 if (controls.target_void_fraction*voidage_rate < total_rate) {
                     actionOnBrokenConstraints(group, Group::InjectionCMode::VREP, reportStepIdx, deferred_logger);
                 }
@@ -1784,6 +1796,9 @@ namespace Opm {
                 sales_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false);
                 sales_rate -= wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true);
 
+                // sum over all nodes
+                sales_rate = comm.sum(sales_rate);
+
                 // add import rate and substract consumption rate for group for gas
                 if (schedule().gConSump(reportStepIdx).has(group.name())) {
                     const auto& gconsump = schedule().gConSump(reportStepIdx).get(group.name(), summaryState);
@@ -1816,6 +1831,10 @@ namespace Opm {
             {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
+
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.oil_target < current_rate  ) {
                     actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::ORAT, reportStepIdx, deferred_logger);
                 }
@@ -1826,6 +1845,9 @@ namespace Opm {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
 
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.water_target < current_rate  ) {
                     actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::WRAT, reportStepIdx, deferred_logger);
                 }
@@ -1834,6 +1856,9 @@ namespace Opm {
             {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
+
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
                 if (controls.gas_target < current_rate  ) {
                     actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::GRAT, reportStepIdx, deferred_logger);
                 }
@@ -1843,6 +1868,10 @@ namespace Opm {
                 double current_rate = 0.0;
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
                 current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
+
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.liquid_target < current_rate  ) {
                     actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::LRAT, reportStepIdx, deferred_logger);
                 }
@@ -1860,6 +1889,9 @@ namespace Opm {
                 current_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true);
                 current_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true);
 
+                // sum over all nodes
+                current_rate = comm.sum(current_rate);
+
                 if (controls.resv_target < current_rate  ) {
                     actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::RESV, reportStepIdx, deferred_logger);
                 }
@@ -1986,6 +2018,10 @@ namespace Opm {
             double gasProductionRate = wellGroupHelpers::sumWellRates(group, schedule, wellState, reportStepIdx, gasPos, /*isInjector*/false);
             double solventProductionRate = wellGroupHelpers::sumSolventRates(group, schedule, wellState, reportStepIdx, /*isInjector*/false);
 
+            const auto& comm = ebosSimulator_.vanguard().grid().comm();
+            solventProductionRate = comm.sum(solventProductionRate);
+            gasProductionRate = comm.sum(gasProductionRate);
+
             double wsolvent = 0.0;
             if (std::abs(gasProductionRate) > 1e-6)
                 wsolvent = solventProductionRate / gasProductionRate;

opm/simulators/wells/WellGroupHelpers.hpp

@@ -295,14 +295,14 @@ namespace Opm {
             wellState.setCurrentProductionGroupReductionRates(group.name(), groupTargetReduction);
     }
 
-
-    inline void updateGuideRateForGroups(const Group& group, const Schedule& schedule, const PhaseUsage& pu, const int reportStepIdx, const double& simTime, const bool isInjector, WellStateFullyImplicitBlackoil& wellState, GuideRate* guideRate, std::vector<double>& pot)
+    template <class Comm>
+    inline void updateGuideRateForGroups(const Group& group, const Schedule& schedule, const PhaseUsage& pu, const int reportStepIdx, const double& simTime, const bool isInjector, WellStateFullyImplicitBlackoil& wellState, const Comm& comm, GuideRate* guideRate, std::vector<double>& pot)
     {
         const int np = pu.num_phases;
         for (const std::string& groupName : group.groups()) {
             std::vector<double> thisPot(np, 0.0);
             const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
-            updateGuideRateForGroups(groupTmp, schedule, pu, reportStepIdx, simTime, isInjector, wellState, guideRate, thisPot);
+            updateGuideRateForGroups(groupTmp, schedule, pu, reportStepIdx, simTime, isInjector, wellState, comm, guideRate, thisPot);
 
             // accumulate group contribution from sub group if FLD
             if (isInjector) {
@@ -354,8 +354,8 @@ namespace Opm {
 
         double oilPot = 0.0;
         if (pu.phase_used[BlackoilPhases::Liquid])
-
             oilPot = pot [ pu.phase_pos[BlackoilPhases::Liquid]];
+
         double gasPot = 0.0;
         if (pu.phase_used[BlackoilPhases::Vapour])
             gasPot = pot [ pu.phase_pos[BlackoilPhases::Vapour]];
@@ -364,6 +364,10 @@ namespace Opm {
         if (pu.phase_used[BlackoilPhases::Aqua])
             waterPot = pot [pu.phase_pos[BlackoilPhases::Aqua]];
 
+        oilPot = comm.sum(oilPot);
+        gasPot = comm.sum(gasPot);
+        waterPot = comm.sum(waterPot);
+
         if (isInjector) {
             wellState.setCurrentGroupInjectionPotentials(group.name(), pot);
         } else {
@@ -371,6 +375,39 @@ namespace Opm {
         }
     }
 
+    template <class Comm>
+    inline void updateGuideRatesForWells(const Schedule& schedule, const PhaseUsage& pu, const int reportStepIdx, const double& simTime, const WellStateFullyImplicitBlackoil& wellState, const Comm& comm, GuideRate* guideRate) {
+
+        const auto& end = wellState.wellMap().end();
+        for (const auto& well : schedule.getWells(reportStepIdx)) {
+            double oilpot = 0.0;
+            double gaspot = 0.0;
+            double waterpot = 0.0;
+
+            const auto& it = wellState.wellMap().find( well.name());
+            if (it != end) {  // the well is found
+
+                int well_index = it->second[0];
+
+                const auto wpot = wellState.wellPotentials().data() + well_index*wellState.numPhases();
+                if (pu.phase_used[BlackoilPhases::Liquid] > 0)
+                    oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]];
+
+                if (pu.phase_used[BlackoilPhases::Vapour] > 0)
+                    gaspot = wpot[pu.phase_pos[BlackoilPhases::Vapour]];
+
+                if (pu.phase_used[BlackoilPhases::Aqua] > 0)
+                    waterpot = wpot[pu.phase_pos[BlackoilPhases::Aqua]];
+            }
+            oilpot = comm.sum(oilpot);
+            gaspot = comm.sum(gaspot);
+            waterpot = comm.sum(waterpot);
+            guideRate->compute(well.name(), reportStepIdx, simTime, oilpot, gaspot, waterpot);
+        }
+
+    }
+
+
     inline void updateVREPForGroups(const Group& group, const Schedule& schedule, const int reportStepIdx, const WellStateFullyImplicitBlackoil& wellStateNupcol, WellStateFullyImplicitBlackoil& wellState, double& resv) {
         for (const std::string& groupName : group.groups()) {
             const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
@@ -415,9 +452,12 @@ namespace Opm {
     inline double wellFractionFromGuideRates(const Well& well, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const GuideRate* guideRate, const Well::GuideRateTarget& wellTarget, const bool isInjector) {
         double groupTotalGuideRate = 0.0;
         const Group& groupTmp = schedule.getGroup(well.groupName(), reportStepIdx);
+        int global_well_index = -1;
         for (const std::string& wellName : groupTmp.wells()) {
             const auto& wellTmp = schedule.getWell(wellName, reportStepIdx);
 
+            global_well_index++;
+
             if (wellTmp.isProducer() && isInjector)
                  continue;
 
@@ -427,19 +467,12 @@ namespace Opm {
             if (wellTmp.getStatus() == Well::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];
-
             // only count wells under group control
             if (isInjector) {
-                if (wellState.currentInjectionControls()[well_index] != Well::InjectorCMode::GRUP)
+                if (!wellState.isInjectionGrup(wellName))
                     continue;
             } else {
-                if (wellState.currentProductionControls()[well_index] !=  Well::ProducerCMode::GRUP)
+                if (!wellState.isProductionGrup(wellName))
                     continue;
             }
 

opm/simulators/wells/WellState.hpp

@@ -63,6 +63,7 @@ namespace Opm
                 const int nw = wells_ecl.size();
                 // const int np = wells->number_of_phases;
                 const int np = pu.num_phases;
+                np_ = np;
                 open_for_output_.assign(nw, true);
                 bhp_.resize(nw, 0.0);
                 thp_.resize(nw, 0.0);
@@ -150,7 +151,7 @@ namespace Opm
         /// The number of phases present.
         int numPhases() const
         {
-            return wellRates().size() / numWells();
+            return np_;
         }
 
 
@@ -222,6 +223,7 @@ namespace Opm
         std::vector<double> wellrates_;
         std::vector<double> perfrates_;
         std::vector<double> perfpress_;
+        int np_;
     protected:
         std::vector<bool>   open_for_output_;
     private:

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -74,11 +74,16 @@ namespace Opm
                   const WellStateFullyImplicitBlackoil* prevState,
                   const PhaseUsage& pu,
                   const std::vector<std::vector<PerforationData>>& well_perf_data,
-                  const SummaryState& summary_state)
+                  const SummaryState& summary_state,
+                  const int globalNumberOfWells)
         {
             // call init on base class
             BaseType :: init(cellPressures, wells_ecl, pu, well_perf_data, summary_state);
 
+            globalIsInjectionGrup_.assign(globalNumberOfWells,0);
+            globalIsProductionGrup_.assign(globalNumberOfWells,0);
+            wellNameToGlobalIdx_.clear();
+
             const int nw = wells_ecl.size();
 
             if( nw == 0 ) return ;
@@ -286,10 +291,11 @@ namespace Opm
                     const size_t numCells,
                     const PhaseUsage& pu,
                     const std::vector<std::vector<PerforationData>>& well_perf_data,
-                    const SummaryState& summary_state)
+                    const SummaryState& summary_state,
+                    const int globalNumWells)
         {
             const std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
-            init(tmp, schedule, wells_ecl, 0, nullptr, pu, well_perf_data, summary_state);
+            init(tmp, schedule, wells_ecl, 0, nullptr, pu, well_perf_data, summary_state, globalNumWells);
 
             if (handle_ms_well) {
                 initWellStateMSWell(wells_ecl, pu, nullptr);
@@ -829,11 +835,82 @@ namespace Opm
                 this->well_reservoir_rates_[np * well_index + p] = 0;
         }
 
+        template<class Comm>
+        void communicateGroupReductionRates(const Comm& comm) {
+            // sum over all nodes
+            for (auto& x : production_group_reduction_rates) {
+                comm.sum(x.second.data(), x.second.size());
+            }
+            for (auto& x : injection_group_reduction_rates) {
+                comm.sum(x.second.data(), x.second.size());
+            }
+        }
+
+        template<class Comm>
+        void communicateReinVrep(const Comm& comm) {
+            // sum over all nodes
+            for (auto& x : injection_group_rein_rates) {
+                comm.sum(x.second.data(), x.second.size());
+            }
+            for (auto& x : injection_group_vrep_rates) {
+                x.second = comm.sum(x.second);
+            }
+        }
+
+        template<class Comm>
+        void updateGlobalIsGrup(const Schedule& schedule, const int reportStepIdx, const Comm& comm) {
+
+            int global_well_index = -1;
+            const auto& end = wellMap().end();
+            for (const auto& well : schedule.getWells(reportStepIdx)) {
+                global_well_index ++;
+                wellNameToGlobalIdx_[well.name()] = global_well_index;
+
+                const auto& it = wellMap().find( well.name());
+                if (it == end)  // the well is not found
+                    continue;
+
+                int well_index = it->second[0];
+
+                if (well.isInjector())
+                    globalIsInjectionGrup_[global_well_index] = (currentInjectionControls()[well_index] == Well::InjectorCMode::GRUP);
+                else
+                    globalIsProductionGrup_[global_well_index] = (currentProductionControls()[well_index] == Well::ProducerCMode::GRUP);
+            }
+            comm.sum(globalIsInjectionGrup_.data(), globalIsInjectionGrup_.size());
+            comm.sum(globalIsProductionGrup_.data(), globalIsProductionGrup_.size());
+        }
+
+        bool isInjectionGrup(const std::string& name) const {
+
+            auto it = wellNameToGlobalIdx_.find(name);
+
+            if (it == wellNameToGlobalIdx_.end())
+                OPM_THROW(std::logic_error, "Could not find global injection group for well" << name);
+
+            return globalIsInjectionGrup_[it->second];
+        }
+
+        bool isProductionGrup(const std::string& name) const {
+
+            auto it = wellNameToGlobalIdx_.find(name);
+
+            if (it == wellNameToGlobalIdx_.end())
+                OPM_THROW(std::logic_error, "Could not find global injection group for well" << name);
+
+            return globalIsProductionGrup_[it->second];
+        }
 
     private:
         std::vector<double> perfphaserates_;
         std::vector<Opm::Well::InjectorCMode> current_injection_controls_;
         std::vector<Well::ProducerCMode> current_production_controls_;
+
+        // size of global number of wells
+        std::vector<int> globalIsInjectionGrup_;
+        std::vector<int> globalIsProductionGrup_;
+        std::map<std::string, int> wellNameToGlobalIdx_;
+
         std::map<std::string, Group::ProductionCMode> current_production_group_controls_;
         std::map<std::string, Group::InjectionCMode> current_injection_group_controls_;
 

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -121,7 +121,7 @@ namespace {
 
         state.init(cpress, setup.sched,
                    setup.sched.getWells(timeStep),
-                   timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st);
+                   timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st, setup.sched.getWells(timeStep).size());
 
         state.initWellStateMSWell(setup.sched.getWells(timeStep),
                                   setup.pu, nullptr);