diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 2f850db88..70530ad69 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -228,7 +228,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 = typename std::remove_reference::type; + 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 ); @@ -247,25 +250,43 @@ 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 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; @@ -285,6 +306,33 @@ namespace Opm WellMapType wellMap_; + using MPIComm = typename Dune::MPIHelper::MPICommunicator; +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7) + using Communication = Dune::Communication; +#else + using Communication = Dune::CollectiveCommunication; +#endif + void gatherVectorsOnRoot(std::vector< data::Connection > from_connections, + std::vector< data::Connection > to_connections, + const Communication& comm) const + { + int size = from_connections.size(); + std::vector sizes; + std::vector 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()); + to_connections.resize(displ.back()); + } + comm.gatherv(from_connections.data(), size, to_connections.data(), + sizes.data(), displ.data(), 0); + } void initSingleWell(const std::vector& cellPressures, const int w, const Well& well, diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 22945c5ed..eb53e673d 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -533,20 +533,16 @@ namespace Opm using rt = data::Rates::opt; std::vector< rt > phs( np ); - std::vector 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 @@ -632,31 +628,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); @@ -668,6 +639,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 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& wells_ecl,