Merge pull request #2954 from blattms/distributed-well-report

2nd part distributed stdwells: Write reports only for wells owned.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -106,7 +106,6 @@ namespace Opm {
 
         // Create cartesian to compressed mapping
         const auto& schedule_wells = schedule().getWellsatEnd();
-        const auto& cartesianSize = Opm::UgGridHelpers::cartDims(grid());
 
         // initialize the additional cell connections introduced by wells.
         for (const auto& well : schedule_wells)
@@ -119,11 +118,8 @@ namespace Opm {
             for ( size_t c=0; c < connectionSet.size(); c++ )
             {
                 const auto& connection = connectionSet.get(c);
-                int i = connection.getI();
-                int j = connection.getJ();
-                int k = connection.getK();
-                int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
-                int compressed_idx = cartesian_to_compressed_.at(cart_grid_idx);
+                int compressed_idx = cartesian_to_compressed_
+                    .at(connection.global_index());
 
                 if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
                     wellCells.push_back(compressed_idx);
@@ -563,6 +559,7 @@ namespace Opm {
         int globalNumWells = 0;
         // Make wells_ecl_ contain only this partition's non-shut wells.
         wells_ecl_ = getLocalNonshutWells(report_step, globalNumWells);
+        local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
 
         this->initializeWellProdIndCalculators();
         initializeWellPerfData();
@@ -572,7 +569,7 @@ 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, globalNumWells); // Resize for restart step
+            well_state_.resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState, globalNumWells); // Resize for restart step
             wellsToState(restartValues.wells, restartValues.grp_nwrk, phaseUsage, handle_ms_well, well_state_);
         }
 
@@ -606,27 +603,26 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     initializeWellPerfData()
     {
-        const auto& grid = ebosSimulator_.vanguard().grid();
-        const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
         well_perf_data_.resize(wells_ecl_.size());
         int well_index = 0;
         for (const auto& well : wells_ecl_) {
             std::size_t completion_index = 0;
             well_perf_data_[well_index].clear();
             well_perf_data_[well_index].reserve(well.getConnections().size());
+            CheckDistributedWellConnections checker(well, *local_parallel_well_info_[well_index]);
+            bool hasFirstPerforation = false;
+            bool firstOpenCompletion = true;
+
             for (const auto& completion : well.getConnections()) {
+                const int active_index =
+                    cartesian_to_compressed_[completion.global_index()];
                 if (completion.state() == Connection::State::OPEN) {
-                    const int i = completion.getI();
-                    const int j = completion.getJ();
-                    const int k = completion.getK();
-                    const int cart_grid_indx = i + cartDims[0] * (j + cartDims[1] * k);
-                    const int active_index = cartesian_to_compressed_[cart_grid_indx];
-                    if (active_index < 0) {
-                        const std::string msg
-                            = ("Cell with i,j,k indices " + std::to_string(i) + " " + std::to_string(j) + " "
-                               + std::to_string(k) + " not found in grid (well = " + well.name() + ").");
-                        OPM_THROW(std::runtime_error, msg);
-                    } else {
+                    if (active_index >= 0) {
+                        if (firstOpenCompletion)
+                        {
+                            hasFirstPerforation = true;
+                        }
+                        checker.connectionFound(completion_index);
                         PerforationData pd;
                         pd.cell_index = active_index;
                         pd.connection_transmissibility_factor = completion.CF();
@@ -634,7 +630,9 @@ namespace Opm {
                         pd.ecl_index = completion_index;
                         well_perf_data_[well_index].push_back(pd);
                     }
+                    firstOpenCompletion = false;
                 } else {
+                    checker.connectionFound(completion_index);
                     if (completion.state() != Connection::State::SHUT) {
                         OPM_THROW(std::runtime_error,
                                   "Completion state: " << Connection::State2String(completion.state()) << " not handled");
@@ -642,6 +640,8 @@ namespace Opm {
                 }
                 ++completion_index;
             }
+            checker.checkAllConnectionsFound();
+            local_parallel_well_info_[well_index]->communicateFirstPerforation(hasFirstPerforation);
             ++well_index;
         }
     }
@@ -686,7 +686,7 @@ namespace Opm {
             }
         }
 
-        well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx,
+        well_state_.init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
                          &previous_well_state_, phase_usage_, well_perf_data_,
                          summaryState, globalNumWells);
     }
@@ -839,12 +839,18 @@ namespace Opm {
         // Use the pvtRegionIdx from the top cell
         const auto& perf_data = this->well_perf_data_[wellID];
 
+        // Cater for case where local part might have no perforations.
+        const int pvtreg = perf_data.empty() ?
+            0 : pvt_region_idx_[perf_data.front().cell_index];
+        const auto& parallel_well_info = *local_parallel_well_info_[wellID];
+        auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
+
         return std::make_unique<WellType>(this->wells_ecl_[wellID],
-                                          *local_parallel_well_info_[wellID],
+                                          parallel_well_info,
                                           time_step,
                                           this->param_,
                                           *this->rateConverter_,
-                                          this->pvt_region_idx_[perf_data.front().cell_index],
+                                          global_pvtreg,
                                           this->numComponents(),
                                           this->numPhases(),
                                           wellID,

opm/simulators/wells/WellGroupHelpers.cpp

@@ -112,6 +112,15 @@ namespace WellGroupHelpers
                 schedule.getGroup(group.parent(), reportStepIdx), schedule, reportStepIdx, factor);
     }
 
+    bool wellIsOwned(std::size_t well_index, [[maybe_unused]] const std::string& wellName,
+                     const WellStateFullyImplicitBlackoil& wellState)
+    {
+        const auto& well_info = wellState.parallelWellInfo(well_index);
+        assert(well_info.name() == wellName);
+
+        return well_info.isOwner();
+    }
+
     double sumWellPhaseRates(const std::vector<double>& rates,
                              const Group& group,
                              const Schedule& schedule,
@@ -135,6 +144,11 @@ namespace WellGroupHelpers
 
             int well_index = it->second[0];
 
+            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            {
+                continue;
+            }
+
             const auto& wellEcl = schedule.getWell(wellName, reportStepIdx);
             // only count producers or injectors
             if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
@@ -195,6 +209,11 @@ namespace WellGroupHelpers
 
             int well_index = it->second[0];
 
+            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            {
+                continue;
+            }
+
             const auto& wellEcl = schedule.getWell(wellName, reportStepIdx);
             // only count producers or injectors
             if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
@@ -303,6 +322,12 @@ namespace WellGroupHelpers
                 continue;
 
             int well_index = it->second[0];
+
+            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            {
+                continue;
+            }
+
             const auto wellrate_index = well_index * wellState.numPhases();
             const double efficiency = wellTmp.getEfficiencyFactor();
             // add contributino from wells not under group control
@@ -385,6 +410,12 @@ namespace WellGroupHelpers
                     continue;
 
                 int well_index = it->second[0];
+
+                if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+                {
+                    continue;
+                }
+
                 const auto wellrate_index = well_index * wellState.numPhases();
                 // add contribution from wells unconditionally
                 for (int phase = 0; phase < np; phase++) {
@@ -431,20 +462,26 @@ namespace WellGroupHelpers
                 double waterpot = 0.0;
 
                 const auto& it = wellState.wellMap().find( well.name());
-                if (it != end) {  // the well is found
+                if (it == end)   // the well is not found
+                    continue;
+                int well_index = it->second[0];
 
-                    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]];
+                if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+                {
+                    continue;
                 }
+
+
+                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]];
+
                 const double wefac = well.getEfficiencyFactor();
                 oilpot = comm.sum(oilpot) * wefac;
                 gaspot = comm.sum(gaspot) * wefac;

opm/simulators/wells/WellInterface_impl.hpp

@@ -1178,7 +1178,6 @@ namespace Opm
                                Opm::DeferredLogger& deferred_logger
                                )
     {
-        const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
         auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
         auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
 
@@ -1194,23 +1193,18 @@ namespace Opm
 
         // COMPDAT handling
         const auto& connectionSet = well_ecl_.getConnections();
+        CheckDistributedWellConnections checker(well_ecl_, parallel_well_info_);
         for (size_t c=0; c<connectionSet.size(); c++) {
             const auto& connection = connectionSet.get(c);
+            const int cell =
+                cartesian_to_compressed[connection.global_index()];
+            if (connection.state() != Connection::State::OPEN || cell >= 0)
+            {
+                checker.connectionFound(c);
+            }
             if (connection.state() == Connection::State::OPEN) {
-                const int i = connection.getI();
-                const int j = connection.getJ();
-                const int k = connection.getK();
 
-                const int* cpgdim = cart_dims;
-                const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
-                const int cell = cartesian_to_compressed[cart_grid_indx];
-
-                if (cell < 0) {
-                    OPM_DEFLOG_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
-                              << k << " not found in grid (well = " << name() << ')', deferred_logger);
-                }
-
-                {
+                if (cell >= 0) {
                     double radius = connection.rw();
                     const std::array<double, 3> cubical =
                         wellhelpers::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
@@ -1242,6 +1236,7 @@ namespace Opm
                 }
             }
         }
+        checker.checkAllConnectionsFound();
     }
 
     template<typename TypeTag>

opm/simulators/wells/WellState.hpp

@@ -25,6 +25,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 #include <opm/simulators/wells/PerforationData.hpp>
+#include <opm/simulators/wells/ParallelWellInfo.hpp>
 
 #include <array>
 #include <cassert>
@@ -51,6 +52,7 @@ namespace Opm
         /// with -1e100.
         void init(const std::vector<double>& cellPressures,
                   const std::vector<Well>& wells_ecl,
+                  const std::vector<ParallelWellInfo*>& parallel_well_info,
                   const PhaseUsage& pu,
                   const std::vector<std::vector<PerforationData>>& well_perf_data,
                   const SummaryState& summary_state)
@@ -59,6 +61,7 @@ namespace Opm
             wellMap_.clear();
 
             well_perf_data_ = well_perf_data;
+            parallel_well_info_ = parallel_well_info;
 
             {
                 // const int nw = wells->number_of_wells;
@@ -76,7 +79,7 @@ namespace Opm
                     const Well& well = wells_ecl[w];
               
                     // Initialize bhp(), thp(), wellRates(), temperature().
-                    initSingleWell(cellPressures, w, well, pu, summary_state);
+                    initSingleWell(cellPressures, w, well, *parallel_well_info[w], pu, summary_state);
 
                     // Setup wellname -> well index mapping.
                     const int num_perf_this_well = well_perf_data[w].size();
@@ -192,6 +195,11 @@ namespace Opm
         const WellMapType& wellMap() const { return wellMap_; }
         WellMapType& wellMap() { return wellMap_; }
 
+        const ParallelWellInfo& parallelWellInfo(std::size_t well_index) const
+        {
+            return *parallel_well_info_[well_index];
+        }
+
         /// The number of wells present.
         int numWells() const
         {
@@ -225,7 +233,10 @@ namespace Opm
                 if (!this->open_for_output_[well_index])
                     continue;
 
-                auto& well = dw[ itr.first ];
+                const auto& pwinfo = *parallel_well_info_[well_index];
+                using WellT = std::remove_reference_t<decltype(dw[ itr.first ])>;
+                WellT dummyWell; // dummy if we are not owner
+                auto& well = pwinfo.isOwner() ? dw[ itr.first ] : dummyWell;
                 well.bhp = this->bhp().at( well_index );
                 well.thp = this->thp().at( well_index );
                 well.temperature = this->temperature().at( well_index );
@@ -244,25 +255,44 @@ namespace Opm
                     well.rates.set( rt::gas, wv[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ] );
                 }
 
-                const auto& pd = this->well_perf_data_[well_index];
-                const int num_perf_well = pd.size();
-                well.connections.resize(num_perf_well);
-
-                for( int i = 0; i < num_perf_well; ++i ) {
-                    const auto active_index = this->well_perf_data_[well_index][i].cell_index;
-                    auto& connection = well.connections[ i ];
-                    connection.index = globalCellIdxMap[active_index];
-                    connection.pressure = this->perfPress()[ itr.second[1] + i ];
-                    connection.reservoir_rate = this->perfRates()[ itr.second[1] + i ];
-                    connection.trans_factor = pd[i].connection_transmissibility_factor;
+                if (pwinfo.communication().size()==1)
+                {
+                    reportConnections(well, pu, itr, globalCellIdxMap);
+                }
+                else
+                {
+                    assert(pwinfo.communication().rank() != 0 || &dummyWell != &well);
+                    // report the local connections
+                    reportConnections(dummyWell, pu, itr, globalCellIdxMap);
+                    // gather them to well on root.
+                    gatherVectorsOnRoot(dummyWell.connections, well.connections,
+                                        pwinfo.communication());
                 }
-                assert(num_perf_well == int(well.connections.size()));
             }
 
             return dw;
 
         }
 
+        virtual void reportConnections(data::Well& well, [[maybe_unused]] const PhaseUsage& pu,
+                                       const WellMapType::value_type& itr,
+                                       const int* globalCellIdxMap) const
+        {
+            const auto well_index = itr.second[ 0 ];
+            const auto& pd = this->well_perf_data_[well_index];
+            const int num_perf_well = pd.size();
+            well.connections.resize(num_perf_well);
+
+            for( int i = 0; i < num_perf_well; ++i ) {
+                const auto active_index = this->well_perf_data_[well_index][i].cell_index;
+                auto& connection = well.connections[ i ];
+                connection.index = globalCellIdxMap[active_index];
+                connection.pressure = this->perfPress()[ itr.second[1] + i ];
+                connection.reservoir_rate = this->perfRates()[ itr.second[1] + i ];
+                connection.trans_factor = pd[i].connection_transmissibility_factor;
+            }
+            assert(num_perf_well == int(well.connections.size()));
+        }
         virtual ~WellState() = default;
         WellState() = default;
         WellState(const WellState& rhs)  = default;
@@ -282,9 +312,37 @@ namespace Opm
 
         WellMapType wellMap_;
 
+        using MPIComm = typename Dune::MPIHelper::MPICommunicator;
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
+        using Communication = Dune::Communication<MPIComm>;
+#else
+        using Communication = Dune::CollectiveCommunication<MPIComm>;
+#endif
+        void gatherVectorsOnRoot(const std::vector< data::Connection >& from_connections,
+                                 std::vector< data::Connection >& to_connections,
+                                 const Communication& comm) const
+        {
+            int size = from_connections.size();
+            std::vector<int> sizes;
+            std::vector<int> displ;
+            if (comm.rank()==0){
+                sizes.resize(comm.size());
+            }
+            comm.gather(&size, sizes.data(), 1, 0);
+
+            if (comm.rank()==0){
+                displ.resize(comm.size()+1, 0);
+                std::transform(displ.begin(), displ.end()-1, sizes.begin(), displ.begin()+1,
+                               std::plus<int>());
+                to_connections.resize(displ.back());
+            }
+            comm.gatherv(from_connections.data(), size, to_connections.data(),
+                         sizes.data(), displ.data(), 0);
+        }
         void initSingleWell(const std::vector<double>& cellPressures,
                             const int w,
                             const Well& well,
+                            const ParallelWellInfo& well_info,
                             const PhaseUsage& pu,
                             const SummaryState& summary_state)
         {
@@ -319,7 +377,9 @@ namespace Opm
                 : (prod_controls.cmode == Well::ProducerCMode::GRUP);
 
             const double inj_surf_rate = well.isInjector() ? inj_controls.surface_rate : 0.0; // To avoid a "maybe-uninitialized" warning.
-
+            const double local_pressure = well_perf_data_[w].empty() ?
+                0 : cellPressures[well_perf_data_[w][0].cell_index];
+            const double global_pressure = well_info.broadcastFirstPerforationValue(local_pressure);
             if (well.getStatus() == Well::Status::STOP) {
                 // Stopped well:
                 // 1. Rates: zero well rates.
@@ -329,8 +389,7 @@ namespace Opm
                 if (is_bhp) {
                     bhp_[w] = bhp_limit;
                 } else {
-                    const int first_cell = well_perf_data_[w][0].cell_index;
-                    bhp_[w] = cellPressures[first_cell];
+                    bhp_[w] = global_pressure;
                 }
             } else if (is_grup) {
                 // Well under group control.
@@ -339,9 +398,8 @@ namespace Opm
                 //    little above or below (depending on if
                 //    the well is an injector or producer)
                 //    pressure in first perforation cell.
-                const int first_cell = well_perf_data_[w][0].cell_index;
                 const double safety_factor = well.isInjector() ? 1.01 : 0.99;
-                bhp_[w] = safety_factor*cellPressures[first_cell];
+                bhp_[w] = safety_factor * global_pressure;
             } else {
                 // Open well, under own control:
                 // 1. Rates: initialize well rates to match
@@ -400,9 +458,8 @@ namespace Opm
                 if (is_bhp) {
                     bhp_[w] = bhp_limit;
                 } else {
-                    const int first_cell = well_perf_data_[w][0].cell_index;
                     const double safety_factor = well.isInjector() ? 1.01 : 0.99;
-                    bhp_[w] = safety_factor*cellPressures[first_cell];
+                    bhp_[w] = safety_factor * global_pressure;
                 }
             }
 
@@ -419,6 +476,7 @@ namespace Opm
 
     protected:
         std::vector<std::vector<PerforationData>> well_perf_data_;
+        std::vector<ParallelWellInfo*> parallel_well_info_;
     };
 
 } // namespace Opm

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -72,6 +72,7 @@ namespace Opm
         void init(const std::vector<double>& cellPressures,
                   const Schedule& schedule,
                   const std::vector<Well>& wells_ecl,
+                  const std::vector<ParallelWellInfo*>& parallel_well_info,
                   const int report_step,
                   const WellStateFullyImplicitBlackoil* prevState,
                   const PhaseUsage& pu,
@@ -80,8 +81,12 @@ namespace Opm
                   const int globalNumberOfWells)
         {
             // call init on base class
-            BaseType :: init(cellPressures, wells_ecl, pu, well_perf_data, summary_state);
+            BaseType :: init(cellPressures, wells_ecl, parallel_well_info, pu, well_perf_data, summary_state);
 
+            for (const auto& winfo: parallel_well_info)
+            {
+                well_rates.insert({winfo->name(), std::make_pair(winfo->isOwner(), std::vector<double>())});
+            }
             globalIsInjectionGrup_.assign(globalNumberOfWells,0);
             globalIsProductionGrup_.assign(globalNumberOfWells,0);
             wellNameToGlobalIdx_.clear();
@@ -319,6 +324,7 @@ namespace Opm
 
 
         void resize(const std::vector<Well>& wells_ecl,
+                    const std::vector<ParallelWellInfo*>& parallel_well_info,
                     const Schedule& schedule,
                     const bool handle_ms_well,
                     const size_t numCells,
@@ -328,7 +334,7 @@ namespace Opm
                     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, globalNumWells);
+            init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, pu, well_perf_data, summary_state, globalNumWells);
 
             if (handle_ms_well) {
                 initWellStateMSWell(wells_ecl, pu, nullptr);
@@ -385,7 +391,7 @@ namespace Opm
         }
 
         void setCurrentWellRates(const std::string& wellName, const std::vector<double>& rates ) {
-            well_rates[wellName] = rates;
+            well_rates[wellName].second = rates;
         }
 
         const std::vector<double>& currentWellRates(const std::string& wellName) const {
@@ -394,7 +400,7 @@ namespace Opm
             if (it == well_rates.end())
                 OPM_THROW(std::logic_error, "Could not find any rates for well  " << wellName);
 
-            return it->second;
+            return it->second.second;
         }
 
         bool hasWellRates(const std::string& wellName) const {
@@ -527,20 +533,16 @@ namespace Opm
 
             using rt = data::Rates::opt;
             std::vector< rt > phs( np );
-            std::vector<rt> pi(np);
             if( pu.phase_used[Water] ) {
                 phs.at( pu.phase_pos[Water] ) = rt::wat;
-                pi .at( pu.phase_pos[Water] ) = rt::productivity_index_water;
             }
 
             if( pu.phase_used[Oil] ) {
                 phs.at( pu.phase_pos[Oil] ) = rt::oil;
-                pi .at( pu.phase_pos[Oil] ) = rt::productivity_index_oil;
             }
 
             if( pu.phase_used[Gas] ) {
                 phs.at( pu.phase_pos[Gas] ) = rt::gas;
-                pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas;
             }
 
             /* this is a reference or example on **how** to convert from
@@ -552,7 +554,8 @@ namespace Opm
 
             for( const auto& wt : this->wellMap() ) {
                 const auto w = wt.second[ 0 ];
-                if (!this->open_for_output_[w])
+                const auto& pwinfo = *parallel_well_info_[w];
+                if (!this->open_for_output_[w] || !pwinfo.isOwner())
                     continue;
 
                 auto& well = res.at( wt.first );
@@ -626,31 +629,6 @@ namespace Opm
                     curr.inj  = this->currentInjectionControls() [w];
                 }
 
-                size_t local_comp_index = 0;
-                for( auto& comp : well.connections) {
-                    const auto connPhaseOffset = np * (wt.second[1] + local_comp_index);
-
-                    const auto rates  = this->perfPhaseRates().begin() + connPhaseOffset;
-                    const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset;
-
-                    for( int i = 0; i < np; ++i ) {
-                        comp.rates.set( phs[ i ], *(rates  + i) );
-                        comp.rates.set( pi [ i ], *(connPI + i) );
-                    }
-                    if ( pu.has_polymer ) {
-                        comp.rates.set( rt::polymer, this->perfRatePolymer()[wt.second[1] + local_comp_index]);
-                    }
-                    if ( pu.has_brine ) {
-                        comp.rates.set( rt::brine, this->perfRateBrine()[wt.second[1] + local_comp_index]);
-                    }
-                    if ( pu.has_solvent || pu.has_zFraction) {
-                        comp.rates.set( rt::solvent, this->perfRateSolvent()[wt.second[1] + local_comp_index]);
-                    }
-
-                    ++local_comp_index;
-                }
-                assert(local_comp_index == this->well_perf_data_[w].size());
-
                 const auto nseg = this->numSegments(w);
                 for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) {
                     const auto seg_no = this->segmentNumber(w, seg_ix);
@@ -662,6 +640,54 @@ namespace Opm
             return res;
         }
 
+        virtual void reportConnections(data::Well& well, const PhaseUsage &pu,
+                                       const WellMapType::value_type& wt,
+                                       const int* globalCellIdxMap) const
+        {
+            using rt = data::Rates::opt;
+            WellState::reportConnections(well, pu, wt, globalCellIdxMap);
+            const int np = pu.num_phases;
+            size_t local_comp_index = 0;
+            std::vector< rt > phs( np );
+            std::vector<rt> pi(np);
+            if( pu.phase_used[Water] ) {
+                phs.at( pu.phase_pos[Water] ) = rt::wat;
+                pi .at( pu.phase_pos[Water] ) = rt::productivity_index_water;
+            }
+
+            if( pu.phase_used[Oil] ) {
+                phs.at( pu.phase_pos[Oil] ) = rt::oil;
+                pi .at( pu.phase_pos[Oil] ) = rt::productivity_index_oil;
+            }
+
+            if( pu.phase_used[Gas] ) {
+                phs.at( pu.phase_pos[Gas] ) = rt::gas;
+                pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas;
+            }
+            for( auto& comp : well.connections) {
+                const auto connPhaseOffset = np * (wt.second[1] + local_comp_index);
+
+                const auto rates  = this->perfPhaseRates().begin() + connPhaseOffset;
+                const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset;
+
+                for( int i = 0; i < np; ++i ) {
+                    comp.rates.set( phs[ i ], *(rates  + i) );
+                    comp.rates.set( pi [ i ], *(connPI + i) );
+                }
+                if ( pu.has_polymer ) {
+                    comp.rates.set( rt::polymer, this->perfRatePolymer()[wt.second[1] + local_comp_index]);
+                }
+                if ( pu.has_brine ) {
+                    comp.rates.set( rt::brine, this->perfRateBrine()[wt.second[1] + local_comp_index]);
+                }
+                if ( pu.has_solvent ) {
+                    comp.rates.set( rt::solvent, this->perfRateSolvent()[wt.second[1] + local_comp_index]);
+                }
+
+                ++local_comp_index;
+            }
+            assert(local_comp_index == this->well_perf_data_[wt.second[0]].size());
+        }
 
         /// init the MS well related.
         void initWellStateMSWell(const std::vector<Well>& wells_ecl,
@@ -1055,11 +1081,27 @@ namespace Opm
             // Create a function that calls some function
             // for all the individual data items to simplify
             // the further code.
-            auto iterateContainer = [](auto& container, auto& func) {
+            auto iterateContainer = [this](auto& container, auto& func) {
                 for (auto& x : container) {
                     func(x.second);
                 }
             };
+            auto iterateRatesContainer = [this](auto& container, auto& func) {
+                for (auto& x : container) {
+                    if (x.second.first)
+                    {
+                        func(x.second.second);
+                    }
+                    else
+                    {
+                        // We might actually store non-zero values for
+                        // distributed wells even if they are not owned.
+                        std::vector<double> dummyRate;
+                        dummyRate.assign(x.second.second.size(), 0);
+                        func(dummyRate);
+                    }
+                }
+            };
 
             auto forAllGroupData = [&](auto& func) {
                 iterateContainer(injection_group_rein_rates, func);
@@ -1067,7 +1109,7 @@ namespace Opm
                 iterateContainer(injection_group_reduction_rates, func);
                 iterateContainer(injection_group_reservoir_rates, func);
                 iterateContainer(production_group_rates, func);
-                iterateContainer(well_rates, func);
+                iterateRatesContainer(well_rates, func);
             };
 
             // Compute the size of the data.
@@ -1217,7 +1259,7 @@ namespace Opm
         std::map<std::string, Group::ProductionCMode> current_production_group_controls_;
         std::map<std::pair<Opm::Phase, std::string>, Group::InjectionCMode> current_injection_group_controls_;
 
-        std::map<std::string, std::vector<double>> well_rates;
+        std::map<std::string, std::pair<bool, std::vector<double>>> well_rates;
         std::map<std::string, std::vector<double>> production_group_rates;
         std::map<std::string, std::vector<double>> production_group_reduction_rates;
         std::map<std::string, std::vector<double>> injection_group_reduction_rates;

tests/MpiFixture.hpp

@@ -0,0 +1,76 @@
+/*
+  Copyright 2020 OPM-OP AS
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef OPM_SIMULATOR_TEST_MPIFIXTURE_HPP
+#define OPM_SIMULATOR_TEST_MPIFIXTURE_HPP
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <boost/test/unit_test.hpp>
+
+class MPIError {
+public:
+  /** @brief Constructor. */
+  MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
+  /** @brief The error string. */
+  std::string errorstring;
+  /** @brief The mpi error code. */
+  int errorcode;
+};
+
+#ifdef HAVE_MPI
+void MPI_err_handler(MPI_Comm *, int *err_code, ...){
+  char *err_string=new char[MPI_MAX_ERROR_STRING];
+  int err_length;
+  MPI_Error_string(*err_code, err_string, &err_length);
+  std::string s(err_string, err_length);
+  std::cerr << "An MPI Error ocurred:"<<std::endl<<s<<std::endl;
+  delete[] err_string;
+  throw MPIError(s, *err_code);
+}
+#endif
+
+struct MPIFixture
+{
+    MPIFixture()
+    {
+#if HAVE_MPI
+        int m_argc = boost::unit_test::framework::master_test_suite().argc;
+        char** m_argv = boost::unit_test::framework::master_test_suite().argv;
+        helper = &Dune::MPIHelper::instance(m_argc, m_argv);
+#ifdef MPI_2
+        MPI_Comm_create_errhandler(MPI_err_handler, &handler);
+        MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler);
+#else
+        MPI_Errhandler_create(MPI_err_handler, &handler);
+        MPI_Errhandler_set(MPI_COMM_WORLD, handler);
+#endif
+#endif
+    }
+    ~MPIFixture()
+    {
+#if HAVE_MPI
+        MPI_Finalize();
+#endif
+    }
+    Dune::MPIHelper* helper;
+#if HAVE_MPI
+    MPI_Errhandler handler;
+#endif
+};
+#endif

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -21,6 +21,7 @@
 
 #define BOOST_TEST_MODULE WellStateFIBOTest
 
+#include "MpiFixture.hpp"
 #include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/parser/eclipse/Python/Python.hpp>
 
@@ -42,6 +43,8 @@
 #include <cstddef>
 #include <string>
 
+BOOST_GLOBAL_FIXTURE(MPIFixture);
+
 struct Setup
 {
     Setup(const std::string& filename)
@@ -114,7 +117,8 @@ struct Setup
 
 namespace {
     Opm::WellStateFullyImplicitBlackoil
-    buildWellState(const Setup& setup, const std::size_t timeStep)
+    buildWellState(const Setup& setup, const std::size_t timeStep,
+                   std::vector<Opm::ParallelWellInfo>& pinfos)
     {
         auto state  = Opm::WellStateFullyImplicitBlackoil{};
 
@@ -122,9 +126,25 @@ namespace {
             std::vector<double>(setup.grid.c_grid()->number_of_cells,
                                 100.0*Opm::unit::barsa);
 
+        auto wells = setup.sched.getWells(timeStep);
+        pinfos.resize(wells.size());
+        std::vector<Opm::ParallelWellInfo*> ppinfos(wells.size());
+        auto pw = pinfos.begin();
+        auto ppw = ppinfos.begin();
+
+        for (const auto& well : wells)
+        {
+            *pw = {well.name()};
+            *ppw = &(*pw);
+            pw->communicateFirstPerforation(true);
+            ++pw;
+            ++ppw;
+        }
+
         state.init(cpress, setup.sched,
-                   setup.sched.getWells(timeStep),
-                   timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st, setup.sched.getWells(timeStep).size());
+                   wells, ppinfos,
+                   timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st,
+                   wells.size());
 
         state.initWellStateMSWell(setup.sched.getWells(timeStep),
                                   setup.pu, nullptr);
@@ -223,7 +243,8 @@ BOOST_AUTO_TEST_CASE(Linearisation)
     const Setup setup{ "msw.data" };
     const auto tstep = std::size_t{0};
 
-    const auto wstate = buildWellState(setup, tstep);
+    std::vector<Opm::ParallelWellInfo> pinfos;
+    const auto wstate = buildWellState(setup, tstep, pinfos);
 
     BOOST_CHECK_EQUAL(wstate.numSegment(), 6 + 1);
 
@@ -244,7 +265,8 @@ BOOST_AUTO_TEST_CASE(Pressure)
     const Setup setup{ "msw.data" };
     const auto tstep = std::size_t{0};
 
-    auto wstate = buildWellState(setup, tstep);
+    std::vector<Opm::ParallelWellInfo> pinfos;
+    auto wstate = buildWellState(setup, tstep, pinfos);
 
     const auto& wells = setup.sched.getWells(tstep);
     const auto prod01_first = wells[0].name() == "PROD01";
@@ -290,7 +312,8 @@ BOOST_AUTO_TEST_CASE(Rates)
     const Setup setup{ "msw.data" };
     const auto tstep = std::size_t{0};
 
-    auto wstate = buildWellState(setup, tstep);
+    std::vector<Opm::ParallelWellInfo> pinfos;
+    auto wstate = buildWellState(setup, tstep, pinfos);
 
     const auto wells = setup.sched.getWells(tstep);
     const auto prod01_first = wells[0].name() == "PROD01";
@@ -361,7 +384,9 @@ BOOST_AUTO_TEST_CASE(STOP_well)
       also for wells in the STOP state.
     */
     const Setup setup{ "wells_manager_data_wellSTOP.data" };
-    auto wstate = buildWellState(setup, 0);
+
+    std::vector<Opm::ParallelWellInfo> pinfos;
+    auto wstate = buildWellState(setup, 0, pinfos);
     for (const auto& p : wstate.perfPress())
         BOOST_CHECK(p > 0);
 }