Merge pull request #2955 from blattms/stdwell-comm-rebase-split-pwell-info

Last part of current distributed standard wells.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -281,6 +281,8 @@ namespace Opm {
             std::vector<WellInterfacePtr > well_container_;
 
             // map from logically cartesian cell indices to compressed ones
+            // cells not in the interior are not mapped. This deactivates
+            // these for distributed wells and make the distribution non-overlapping.
             std::vector<int> cartesian_to_compressed_;
 
             std::vector<bool> is_cell_perforated_;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -464,6 +464,11 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::
     endReportStep() {
+        // Clear the communication data structures for above values.
+        for(auto&& pinfo : local_parallel_well_info_)
+        {
+            pinfo->clearCommunicateAbove();
+        }
     }
 
     // called at the end of a report step
@@ -614,12 +619,15 @@ namespace Opm {
         well_perf_data_.resize(wells_ecl_.size());
         int well_index = 0;
         for (const auto& well : wells_ecl_) {
-            std::size_t completion_index = 0;
+            int completion_index = 0;
+            int completion_index_above = -1; // -1 marks no above perf available
             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;
+            auto& parallelWellInfo = *local_parallel_well_info_[well_index];
+            parallelWellInfo.beginReset();
 
             for (const auto& completion : well.getConnections()) {
                 const int active_index =
@@ -637,6 +645,8 @@ namespace Opm {
                         pd.satnum_id = completion.satTableId();
                         pd.ecl_index = completion_index;
                         well_perf_data_[well_index].push_back(pd);
+                        parallelWellInfo.pushBackEclIndex(completion_index_above,
+                                                          completion_index);
                     }
                     firstOpenCompletion = false;
                 } else {
@@ -646,10 +656,14 @@ namespace Opm {
                                   "Completion state: " << Connection::State2String(completion.state()) << " not handled");
                     }
                 }
+                // Note: we rely on the connections being filtered! I.e. there are only connections
+                // to active cells in the global grid.
                 ++completion_index;
+                ++completion_index_above;
             }
+            parallelWellInfo.endReset();
             checker.checkAllConnectionsFound();
-            local_parallel_well_info_[well_index]->communicateFirstPerforation(hasFirstPerforation);
+            parallelWellInfo.communicateFirstPerforation(hasFirstPerforation);
             ++well_index;
         }
     }
@@ -1605,8 +1619,15 @@ namespace Opm {
     {
         cartesian_to_compressed_.resize(number_of_cartesian_cells, -1);
         if (global_cell) {
+            auto elemIt = ebosSimulator_.gridView().template begin</*codim=*/ 0>();
             for (unsigned i = 0; i < local_num_cells_; ++i) {
-                cartesian_to_compressed_[global_cell[i]] = i;
+                // Skip perforations in the overlap/ghost for distributed wells.
+                if (elemIt->partitionType() == Dune::InteriorEntity)
+                {
+                    assert(ebosSimulator_.gridView().indexSet().index(*elemIt) == static_cast<int>(i));
+                    cartesian_to_compressed_[global_cell[i]] = i;
+                }
+                ++elemIt;
             }
         }
         else {
@@ -2485,7 +2506,7 @@ namespace Opm {
         // What is the proper approach?
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         const int fipnum = 0;
-        int pvtreg = well_perf_data_.empty()
+        int pvtreg = well_perf_data_.empty() || well_perf_data_[0].empty()
             ? pvt_region_idx_[0]
             : pvt_region_idx_[well_perf_data_[0][0].cell_index];
 

opm/simulators/wells/ParallelWellInfo.hpp

@@ -136,8 +136,17 @@ public:
     template<class T>
     T broadcastFirstPerforationValue(const T& t) const
     {
-        T res=t;
+        T res = t;
+#ifndef NDEBUG
+        assert(rankWithFirstPerf_ >= 0 && rankWithFirstPerf_ < comm_->size());
+        // At least on some OpenMPI version this might broadcast might interfere
+        // with other communication if there are bugs
+        comm_->barrier();
+#endif
         comm_->broadcast(&res, 1, rankWithFirstPerf_);
+#ifndef NDEBUG
+        comm_->barrier();
+#endif
         return res;
     }
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -1880,6 +1880,10 @@ namespace Opm
         }
 
         // Compute the average pressure in each well block
+        auto p_above =  this->parallel_well_info_.communicateAboveValues(well_state.bhp()[w],
+                                                                          well_state.perfPress().data()+first_perf_,
+                                                                          nperf);
+
         for (int perf = 0; perf < nperf; ++perf) {
             const int cell_idx = well_cells_[perf];
             const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
@@ -1887,8 +1891,7 @@ namespace Opm
 
             // TODO: this is another place to show why WellState need to be a vector of WellState.
             // TODO: to check why should be perf - 1
-            const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1];
+            const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above[perf])/2;
-            const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
             const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
             const double saltConcentration = fs.saltConcentration().value();
 
@@ -2126,10 +2129,10 @@ namespace Opm
 
         const int nperf = number_of_perforations_;
         perf_pressure_diffs_.resize(nperf, 0.0);
+        auto z_above = this->parallel_well_info_.communicateAboveValues(ref_depth_, perf_depth_);
 
         for (int perf = 0; perf < nperf; ++perf) {
-            const double z_above = perf == 0 ? ref_depth_ : perf_depth_[perf - 1];
+            const double dz = perf_depth_[perf] - z_above[perf];
-            const double dz = perf_depth_[perf] - z_above;
             perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * gravity_;
         }
 

opm/simulators/wells/WellGroupHelpers.cpp

@@ -113,15 +113,6 @@ 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,
@@ -145,7 +136,7 @@ namespace WellGroupHelpers
 
             int well_index = it->second[0];
 
-            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
             {
                 continue;
             }
@@ -210,7 +201,7 @@ namespace WellGroupHelpers
 
             int well_index = it->second[0];
 
-            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
             {
                 continue;
             }
@@ -324,7 +315,7 @@ namespace WellGroupHelpers
 
             int well_index = it->second[0];
 
-            if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
+            if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
             {
                 continue;
             }
@@ -463,27 +454,22 @@ namespace WellGroupHelpers
                 double waterpot = 0.0;
 
                 const auto& it = wellState.wellMap().find( well.name());
-                if (it == end)   // the well is not found
+                if (it != end && wellState.wellIsOwned(it->second[0], well.name()))
-                    continue;
-                int well_index = it->second[0];
-
-                if (! wellIsOwned(well_index, wellName, wellState) ) // Only sum once
                 {
-                    continue;
+                    // the well is found and owned
+
+                    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();
                 }
-
-
-                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;
                 waterpot = comm.sum(waterpot) * wefac;

opm/simulators/wells/WellGroupHelpers.hpp

@@ -160,6 +160,12 @@ namespace WellGroupHelpers
                 continue;
 
             int well_index = it->second[0];
+
+            if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
+            {
+                continue;
+            }
+
             const auto wellrate_index = well_index * wellState.numPhases();
             // add contribution from wells unconditionally
             for (int phase = 0; phase < np; phase++) {
@@ -207,8 +213,9 @@ namespace WellGroupHelpers
             double waterpot = 0.0;
 
             const auto& it = wellState.wellMap().find(well.name());
-            if (it != end) { // the well is found
+            if (it != end && wellState.wellIsOwned(it->second[0], well.name()))
-
+            {
+                // the well is found and owned
                 int well_index = it->second[0];
 
                 const auto wpot = wellState.wellPotentials().data() + well_index * wellState.numPhases();

opm/simulators/wells/WellInterface_impl.hpp

@@ -859,6 +859,7 @@ namespace Opm
                 }
             } // end of for (const int c : conns)
 
+            parallel_well_info_.communication().sum(completion_rates.data(), completion_rates.size());
             const double ratio_completion = ratioFunc(completion_rates, phaseUsage());
 
             if (ratio_completion > max_ratio_completion) {

opm/simulators/wells/WellState.hpp

@@ -200,6 +200,13 @@ namespace Opm
             return *parallel_well_info_[well_index];
         }
 
+        bool wellIsOwned(std::size_t well_index, [[maybe_unused]] const std::string& wellName) const
+        {
+            const auto& well_info = parallelWellInfo(well_index);
+            assert(well_info.name() == wellName);
+
+            return well_info.isOwner();
+        }
         /// The number of wells present.
         int numWells() const
         {
@@ -332,8 +339,7 @@ namespace Opm
 
             if (comm.rank()==0){
                 displ.resize(comm.size()+1, 0);
-                std::transform(displ.begin(), displ.end()-1, sizes.begin(), displ.begin()+1,
+                std::partial_sum(sizes.begin(), sizes.end(), displ.begin()+1);
-                               std::plus<int>());
                 to_connections.resize(displ.back());
             }
             comm.gatherv(from_connections.data(), size, to_connections.data(),
@@ -359,7 +365,7 @@ namespace Opm
                 temperature_[w] = well.injectionControls(summary_state).temperature;
             }
 
-            const int num_perf_this_well = well_perf_data_[w].size();
+            const int num_perf_this_well = well_info.communication().sum(well_perf_data_[w].size());
             if ( num_perf_this_well == 0 ) {
                 // No perforations of the well. Initialize to zero.
                 bhp_[w] = 0.;