Correct pressure and pvt region index for distributed wells.

To get this we need to determine which process has the
very first open well connection and use that to broadcast
pressure and pvt region index.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -573,7 +573,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_);
         }
 
@@ -616,6 +616,8 @@ namespace Opm {
             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 i = completion.getI();
@@ -625,6 +627,10 @@ namespace Opm {
                 const int active_index = cartesian_to_compressed_[cart_grid_indx];
                 if (completion.state() == Connection::State::OPEN) {
                     if (active_index >= 0) {
+                        if (firstOpenCompletion)
+                        {
+                            hasFirstPerforation = true;
+                        }
                         checker.connectionFound(completion_index);
                         PerforationData pd;
                         pd.cell_index = active_index;
@@ -633,6 +639,7 @@ 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) {
@@ -643,6 +650,7 @@ namespace Opm {
                 ++completion_index;
             }
             checker.checkAllConnectionsFound();
+            local_parallel_well_info_[well_index]->communicateFirstPerforation(hasFirstPerforation);
             ++well_index;
         }
     }
@@ -687,7 +695,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);
     }
@@ -840,12 +848,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.size() ?
+            pvt_region_idx_[perf_data.front().cell_index] : 0;
+        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/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)
@@ -76,7 +78,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();
@@ -285,6 +287,7 @@ namespace Opm
         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 +322,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].size() ?
+                cellPressures[well_perf_data_[w][0].cell_index] : 0;
+            const double global_pressure = well_info.broadcastFirstPerforationValue(local_pressure);
             if (well.getStatus() == Well::Status::STOP) {
                 // Stopped well:
                 // 1. Rates: zero well rates.
@@ -329,8 +334,7 @@ namespace Opm
                 if (is_bhp) {
                     bhp_[w] = bhp_limit;
                 } else {
-                    const int first_cell = well_perf_data_[w][0].cell_index;
+                    bhp_[w] = global_pressure;
-                    bhp_[w] = cellPressures[first_cell];
                 }
             } else if (is_grup) {
                 // Well under group control.
@@ -339,9 +343,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 +403,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;
                 }
             }
 

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,7 +81,7 @@ 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);
 
             globalIsInjectionGrup_.assign(globalNumberOfWells,0);
             globalIsProductionGrup_.assign(globalNumberOfWells,0);
@@ -319,6 +320,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 +330,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);

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -122,8 +122,22 @@ namespace {
             std::vector<double>(setup.grid.c_grid()->number_of_cells,
                                 100.0*Opm::unit::barsa);
 
+        auto wells = setup.sched.getWells(timeStep);
+        std::vector<Opm::ParallelWellInfo> pinfos(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;
+            ++ppw;
+        }
+
         state.init(cpress, setup.sched,
-                   setup.sched.getWells(timeStep),
+                   wells, ppinfos,
                    timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st, setup.sched.getWells(timeStep).size());
 
         state.initWellStateMSWell(setup.sched.getWells(timeStep),