diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 534f187c2..0fd483f00 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -265,7 +265,7 @@ initFromRestartFile(const RestartValue& restartValues, const int report_step = std::max(eclState_.getInitConfig().getRestartStep() - 1, 0); // wells_ecl_ should only contain wells on this processor. wells_ecl_ = getLocalWells(report_step); - local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); + this->local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); this->initializeWellProdIndCalculators(); initializeWellPerfData(); @@ -273,7 +273,7 @@ initFromRestartFile(const RestartValue& restartValues, const int nw = wells_ecl_.size(); if (nw > 0) { handle_ms_well &= anyMSWellOpenLocal(); - this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState_); // Resize for restart step + this->wellState().resize(wells_ecl_, this->local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState_); // Resize for restart step loadRestartData(restartValues.wells, restartValues.grp_nwrk, phase_usage_, handle_ms_well, this->wellState()); } @@ -297,11 +297,11 @@ getLocalWells(const int timeStepIdx) const return w; } -std::vector +std::vector> BlackoilWellModelGeneric:: createLocalParallelWellInfo(const std::vector& wells) { - std::vector local_parallel_well_info; + std::vector> local_parallel_well_info; local_parallel_well_info.reserve(wells.size()); for (const auto& well : wells) { @@ -311,7 +311,7 @@ createLocalParallelWellInfo(const std::vector& wells) wellPair); assert(pwell != parallel_well_info_.end() && *pwell == wellPair); - local_parallel_well_info.push_back(&(*pwell)); + local_parallel_well_info.push_back(std::ref(*pwell)); } return local_parallel_well_info; } @@ -339,10 +339,10 @@ initializeWellPerfData() int completion_index_above = ParallelWellInfo::INVALID_ECL_INDEX; well_perf_data_[well_index].clear(); well_perf_data_[well_index].reserve(well.getConnections().size()); - CheckDistributedWellConnections checker(well, *local_parallel_well_info_[well_index]); + CheckDistributedWellConnections checker(well, local_parallel_well_info_[well_index].get()); bool hasFirstPerforation = false; bool firstOpenCompletion = true; - auto& parallelWellInfo = *local_parallel_well_info_[well_index]; + auto& parallelWellInfo = this->local_parallel_well_info_[well_index].get(); parallelWellInfo.beginReset(); for (const auto& completion : well.getConnections()) { diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 114662f76..11ed887b7 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -241,7 +241,7 @@ protected: /// \brief Create the parallel well information /// \param localWells The local wells from ECL schedule - std::vector createLocalParallelWellInfo(const std::vector& wells); + std::vector> createLocalParallelWellInfo(const std::vector& wells); void initializeWellProdIndCalculators(); void initializeWellPerfData(); @@ -396,7 +396,7 @@ protected: std::vector local_shut_wells_{}; std::vector parallel_well_info_; - std::vector local_parallel_well_info_; + std::vector> local_parallel_well_info_; std::vector prod_index_calc_; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 7bf2ea411..66eee66d6 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -183,7 +183,7 @@ namespace Opm { const auto& summaryState = ebosSimulator_.vanguard().summaryState(); // Make wells_ecl_ contain only this partition's wells. wells_ecl_ = getLocalWells(timeStepIdx); - local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); + this->local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); // at least initializeWellState might be throw // exception in opm-material (UniformTabulated2DFunction.hpp) @@ -416,9 +416,9 @@ namespace Opm { endReportStep() { // Clear the communication data structures for above values. - for (auto&& pinfo : local_parallel_well_info_) + for (auto&& pinfo : this->local_parallel_well_info_) { - pinfo->clear(); + pinfo.get().clear(); } } @@ -736,7 +736,7 @@ namespace Opm { const auto pvtreg = perf_data.empty() ? 0 : pvt_region_idx_[perf_data.front().cell_index]; - const auto& parallel_well_info = *local_parallel_well_info_[wellID]; + const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get(); const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg); return std::make_unique(this->wells_ecl_[wellID], @@ -1685,7 +1685,7 @@ namespace Opm { double weighted_temperature = 0.0; double total_weight = 0.0; - auto& well_info = *local_parallel_well_info_[wellID]; + auto& well_info = local_parallel_well_info_[wellID].get(); const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size()); auto& ws = this->wellState().well(wellID); auto& perf_data = ws.perf_data; diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp index 842abd88e..c3d69d911 100644 --- a/opm/simulators/wells/SingleWellState.cpp +++ b/opm/simulators/wells/SingleWellState.cpp @@ -21,8 +21,9 @@ namespace Opm { -SingleWellState::SingleWellState(bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp) - : producer(is_producer) +SingleWellState::SingleWellState(const ParallelWellInfo& pinfo, bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp) + : parallel_info(pinfo) + , producer(is_producer) , temperature(temp) , well_potentials(num_phases) , productivity_index(num_phases) @@ -66,6 +67,25 @@ void SingleWellState::open() { this->status = Well::Status::OPEN; } + +double SingleWellState::sum_connection_rates(const std::vector connection_rates) const { + return this->parallel_info.get().sumPerfValues(connection_rates.begin(), connection_rates.end()); +} + +double SingleWellState::sum_brine_rates() const { + return this->sum_connection_rates(this->perf_data.brine_rates); +} + + +double SingleWellState::sum_polymer_rates() const { + return this->sum_connection_rates(this->perf_data.polymer_rates); +} + + +double SingleWellState::sum_solvent_rates() const { + return this->sum_connection_rates(this->perf_data.solvent_rates); +} + } diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index 0b09eb120..aa0b079a2 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -20,18 +20,22 @@ #ifndef OPM_SINGLE_WELL_STATE_HEADER_INCLUDED #define OPM_SINGLE_WELL_STATE_HEADER_INCLUDED +#include #include #include #include #include #include +#include namespace Opm { class SingleWellState { public: - SingleWellState(bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp); + SingleWellState(const ParallelWellInfo& pinfo, bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp); + + std::reference_wrapper parallel_info; Well::Status status{Well::Status::OPEN}; bool producer; @@ -55,6 +59,15 @@ public: void shut(); void stop(); void open(); + + // The sum_xxx_rates() functions sum over all connection rates of pertinent + // types. In the case of distributed wells this involves an MPI + // communication. + double sum_solvent_rates() const; + double sum_polymer_rates() const; + double sum_brine_rates() const; +private: + double sum_connection_rates(const std::vector connection_rates) const; }; diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 5a05d08b7..ba081200d 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -291,10 +291,10 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log } if constexpr (has_gfrac_variable) { primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]] - - (Indices::enableSolvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ; + - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ; } if constexpr (Indices::enableSolvent) { - primary_variables_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; + primary_variables_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * ws.sum_solvent_rates() / total_well_rate ; } } else { // total_well_rate == 0 if (baseif_.isInjector()) { diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 6a4a1b356..184a89bcf 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1210,7 +1210,6 @@ namespace Opm const PhaseUsage& pu = phaseUsage(); b_perf.resize(nperf * this->num_components_); surf_dens_perf.resize(nperf * this->num_components_); - const int w = this->index_of_well_; const auto& ws = well_state.well(this->index_of_well_); const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); @@ -1234,8 +1233,6 @@ namespace Opm const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); - // 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_avg = (perf_press[perf] + p_above[perf])/2; const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); const double saltConcentration = fs.saltConcentration().value(); @@ -1254,7 +1251,7 @@ namespace Opm const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); if (oilrate > 0) { - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0); + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); double rv = 0.0; if (gasrate > 0) { rv = oilrate / gasrate; @@ -1277,7 +1274,7 @@ namespace Opm const int oilpos = oilCompIdx + perf * this->num_components_; if (gasPresent) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0); + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); if (gasrate > 0) { const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); double rs = 0.0; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 2f0266831..5d584b733 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -240,11 +240,12 @@ namespace WellGroupHelpers if (wellEcl.getStatus() == Well::Status::SHUT) continue; + const auto& ws = wellState.well(well_index.value()); double factor = wellEcl.getEfficiencyFactor(); if (injector) - rate += factor * wellState.solventWellRate(well_index.value()); + rate += factor * ws.sum_solvent_rates(); else - rate -= factor * wellState.solventWellRate(well_index.value()); + rate -= factor * ws.sum_solvent_rates(); } return rate; } diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 6fa5b3bae..f5854c8e9 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -33,13 +33,12 @@ namespace Opm { void WellState::base_init(const std::vector& cellPressures, - const std::vector& wells_ecl, - const std::vector& parallel_well_info, - const std::vector>& well_perf_data, - const SummaryState& summary_state) + const std::vector& wells_ecl, + const std::vector>& parallel_well_info, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { // clear old name mapping - this->parallel_well_info_.clear(); this->wells_.clear(); { // const int nw = wells->number_of_wells; @@ -61,7 +60,7 @@ void WellState::base_init(const std::vector& cellPressures, void WellState::initSingleWell(const std::vector& cellPressures, const Well& well, const std::vector& well_perf_data, - const ParallelWellInfo* well_info, + const ParallelWellInfo& well_info, const SummaryState& summary_state) { assert(well.isInjector() || well.isProducer()); @@ -72,8 +71,7 @@ void WellState::initSingleWell(const std::vector& cellPressures, const int np = pu.num_phases; double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56; - this->parallel_well_info_.add(well.name(), well_info); - auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), well_perf_data.size(), static_cast(np), temp}); + auto& ws = this->wells_.add(well.name(), SingleWellState{well_info, well.isProducer(), well_perf_data.size(), static_cast(np), temp}); if ( ws.perf_data.empty()) return; @@ -91,7 +89,7 @@ void WellState::initSingleWell(const std::vector& cellPressures, const double local_pressure = well_perf_data.empty() ? 0 : cellPressures[well_perf_data[0].cell_index]; - const double global_pressure = well_info->broadcastFirstPerforationValue(local_pressure); + const double global_pressure = well_info.broadcastFirstPerforationValue(local_pressure); if (well.getStatus() == Well::Status::OPEN) { ws.status = Well::Status::OPEN; @@ -194,13 +192,13 @@ void WellState::initSingleWell(const std::vector& cellPressures, void WellState::init(const std::vector& cellPressures, - const Schedule& schedule, - const std::vector& wells_ecl, - const std::vector& parallel_well_info, - const int report_step, - const WellState* prevState, - const std::vector>& well_perf_data, - const SummaryState& summary_state) + const Schedule& schedule, + const std::vector& wells_ecl, + const std::vector>& parallel_well_info, + const int report_step, + const WellState* prevState, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { // call init on base class this->base_init(cellPressures, wells_ecl, parallel_well_info, well_perf_data, summary_state); @@ -211,7 +209,7 @@ void WellState::init(const std::vector& cellPressures, } for (const auto& winfo: parallel_well_info) { - well_rates[winfo->name()].first = winfo->isOwner(); + well_rates[winfo.get().name()].first = winfo.get().isOwner(); } const int nw = wells_ecl.size(); @@ -371,12 +369,12 @@ void WellState::init(const std::vector& cellPressures, } void WellState::resize(const std::vector& wells_ecl, - const std::vector& parallel_well_info, - const Schedule& schedule, - const bool handle_ms_well, - const size_t numCells, - const std::vector>& well_perf_data, - const SummaryState& summary_state) + const std::vector>& parallel_well_info, + const Schedule& schedule, + const bool handle_ms_well, + const size_t numCells, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { const std::vector tmp(numCells, 0.0); // <- UGLY HACK to pass the size init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, well_perf_data, summary_state); @@ -471,15 +469,15 @@ WellState::report(const int* globalCellIdxMap, } if (pu.has_solvent || pu.has_zFraction) { - well.rates.set(rt::solvent, solventWellRate(well_index)); + well.rates.set(rt::solvent, ws.sum_solvent_rates()); } if (pu.has_polymer) { - well.rates.set(rt::polymer, polymerWellRate(well_index)); + well.rates.set(rt::polymer, ws.sum_polymer_rates()); } if (pu.has_brine) { - well.rates.set(rt::brine, brineWellRate(well_index)); + well.rates.set(rt::brine, ws.sum_brine_rates()); } if (ws.producer) { @@ -500,7 +498,7 @@ WellState::report(const int* globalCellIdxMap, curr.inj = ws.injection_cmode; } - const auto& pwinfo = *this->parallel_well_info_[well_index]; + const auto& pwinfo = ws.parallel_info.get(); if (pwinfo.communication().size()==1) { reportConnections(well.connections, pu, well_index, globalCellIdxMap); @@ -737,31 +735,6 @@ WellState::calculateSegmentRates(const std::vector>& segment_in } } -double WellState::solventWellRate(const int w) const -{ - auto& ws = this->well(w); - const auto& perf_data = ws.perf_data; - const auto& perf_rates_solvent = perf_data.solvent_rates; - return parallel_well_info_[w]->sumPerfValues(perf_rates_solvent.begin(), perf_rates_solvent.end()); -} - -double WellState::polymerWellRate(const int w) const -{ - auto& ws = this->well(w); - const auto& perf_data = ws.perf_data; - const auto& perf_rates_polymer = perf_data.polymer_rates; - return parallel_well_info_[w]->sumPerfValues(perf_rates_polymer.begin(), perf_rates_polymer.end()); -} - -double WellState::brineWellRate(const int w) const -{ - auto& ws = this->well(w); - const auto& perf_data = ws.perf_data; - const auto& perf_rates_brine = perf_data.brine_rates; - return parallel_well_info_[w]->sumPerfValues(perf_rates_brine.begin(), perf_rates_brine.end()); -} - - void WellState::stopWell(int well_index) { auto& ws = this->well(well_index); @@ -908,7 +881,7 @@ WellState::reportSegmentResults(const PhaseUsage& pu, bool WellState::wellIsOwned(std::size_t well_index, [[maybe_unused]] const std::string& wellName) const { - const auto& well_info = parallelWellInfo(well_index); + const auto& well_info = this->parallelWellInfo(well_index); assert(well_info.name() == wellName); return well_info.isOwner(); @@ -986,7 +959,8 @@ void WellState::resetConnectionTransFactors(const int well_index, const ParallelWellInfo& WellState::parallelWellInfo(std::size_t well_index) const { - return *parallel_well_info_[well_index]; + const auto& ws = this->well(well_index); + return ws.parallel_info; } template void WellState::updateGlobalIsGrup(const ParallelWellInfo::Communication& comm); diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 99ea37e01..6290443ff 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -86,14 +86,14 @@ public: void init(const std::vector& cellPressures, const Schedule& schedule, const std::vector& wells_ecl, - const std::vector& parallel_well_info, + const std::vector>& parallel_well_info, const int report_step, const WellState* prevState, const std::vector>& well_perf_data, const SummaryState& summary_state); void resize(const std::vector& wells_ecl, - const std::vector& parallel_well_info, + const std::vector>& parallel_well_info, const Schedule& schedule, const bool handle_ms_well, const size_t numCells, @@ -132,15 +132,6 @@ public: static void calculateSegmentRates(const std::vector>& segment_inlets, const std::vector>&segment_perforations, const std::vector& perforation_rates, const int np, const int segment, std::vector& segment_rates); - /// One rate pr well - double solventWellRate(const int w) const; - - /// One rate pr well - double polymerWellRate(const int w) const; - - /// One rate pr well - double brineWellRate(const int w) const; - template void communicateGroupRates(const Comm& comm); @@ -268,7 +259,6 @@ private: PhaseUsage phase_usage_; WellContainer wells_; - WellContainer parallel_well_info_; // The well_rates variable is defined for all wells on all processors. The // bool in the value pair is whether the current process owns the well or // not. @@ -303,14 +293,14 @@ private: /// with -1e100. void base_init(const std::vector& cellPressures, const std::vector& wells_ecl, - const std::vector& parallel_well_info, + const std::vector>& parallel_well_info, const std::vector>& well_perf_data, const SummaryState& summary_state); void initSingleWell(const std::vector& cellPressures, const Well& well, const std::vector& well_perf_data, - const ParallelWellInfo* well_info, + const ParallelWellInfo& well_info, const SummaryState& summary_state); diff --git a/tests/test_wellstate.cpp b/tests/test_wellstate.cpp index 33150498a..d59c4c2c2 100644 --- a/tests/test_wellstate.cpp +++ b/tests/test_wellstate.cpp @@ -18,6 +18,7 @@ */ #include +#include #define BOOST_TEST_MODULE WellStateFIBOTest @@ -138,17 +139,15 @@ namespace { auto wells = setup.sched.getWells(timeStep); pinfos.resize(wells.size()); - std::vector ppinfos(wells.size()); + std::vector> ppinfos; auto pw = pinfos.begin(); - auto ppw = ppinfos.begin(); for (const auto& well : wells) { *pw = {well.name()}; - *ppw = &(*pw); + ppinfos.push_back(std::ref(*pw)); pw->communicateFirstPerforation(true); ++pw; - ++ppw; } state.init(cpress, setup.sched, @@ -569,9 +568,10 @@ BOOST_AUTO_TEST_CASE(TESTPerfData) { BOOST_AUTO_TEST_CASE(TestSingleWellState) { - Opm::SingleWellState ws1(true, 10, 3, 1); - Opm::SingleWellState ws2(true, 10, 3, 2); - Opm::SingleWellState ws3(false, 10, 3, 3); + Opm::ParallelWellInfo pinfo; + Opm::SingleWellState ws1(pinfo, true, 10, 3, 1); + Opm::SingleWellState ws2(pinfo, true, 10, 3, 2); + Opm::SingleWellState ws3(pinfo, false, 10, 3, 3); ws1.bhp = 100; ws1.thp = 200;