diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 21afaac03..d87d43037 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1619,7 +1619,7 @@ namespace Opm { // potentials is resized and set to zero in the beginning of well->ComputeWellPotentials // and updated only if sucessfull. i.e. the potentials are zero for exceptions for (int p = 0; p < np; ++p) { - this->wellState().wellPotentials()[well->indexOfWell() * np + p] = std::abs(potentials[p]); + this->wellState().wellPotentials(well->indexOfWell())[p] = std::abs(potentials[p]); } } } @@ -2787,8 +2787,7 @@ namespace Opm { wellPI(const int well_index) const { const auto& pu = this->phase_usage_; - const auto np = this->numPhases(); - const auto* pi = &this->wellState().productivityIndex()[np*well_index + 0]; + const auto& pi = this->wellState().productivityIndex(well_index); const auto preferred = this->wells_ecl_[well_index].getPreferredPhase(); switch (preferred) { // Should really have LIQUID = OIL + WATER here too... diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ec60ac383..7569907cb 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -632,7 +632,7 @@ namespace Opm const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0; for (int phase = 0; phase < np; ++phase){ well_state_copy.wellRates(well_copy.index_of_well_)[phase] - = sign * well_state_copy.wellPotentials()[well_copy.index_of_well_*np + phase]; + = sign * well_state_copy.wellPotentials(well_copy.index_of_well_)[phase]; } well_copy.scaleSegmentRatesWithWellRates(well_state_copy); @@ -992,8 +992,8 @@ namespace Opm std::transform(src, src + np, dest, dest, std::plus<>{}); }; - auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0]; - auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0]; + auto* wellPI = well_state.productivityIndex(this->index_of_well_).data(); + auto* connPI = well_state.connectionProductivityIndex(this->index_of_well_).data(); setToZero(wellPI); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 047467d26..511f9fe30 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2042,8 +2042,8 @@ namespace Opm std::transform(src, src + np, dest, dest, std::plus<>{}); }; - auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0]; - auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0]; + auto* wellPI = well_state.productivityIndex(this->index_of_well_).data(); + auto* connPI = well_state.connectionProductivityIndex(this->index_of_well_).data(); setToZero(wellPI); diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index e3e83aae0..b2eac79ca 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -1279,10 +1279,9 @@ namespace WellGroupHelpers continue; } - const auto wellrate_index = well_index * wellState.numPhases(); // add contribution from wells unconditionally for (int phase = 0; phase < np; phase++) { - pot[phase] += wefac * wellState.wellPotentials()[wellrate_index + phase]; + pot[phase] += wefac * wellState.wellPotentials(well_index)[phase]; } } @@ -1332,7 +1331,7 @@ namespace WellGroupHelpers // the well is found and owned int well_index = it->second[0]; - const auto wpot = wellState.wellPotentials().data() + well_index * wellState.numPhases(); + const auto wpot = wellState.wellPotentials(well_index); if (pu.phase_used[BlackoilPhases::Liquid] > 0) oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]]; diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 3fd14321f..26127ba4d 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -314,7 +314,7 @@ protected: double scalingFactor(const int comp_idx) const; - std::vector initialWellRateFractions(const Simulator& ebosSimulator, const std::vector& potentials) const; + std::vector initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const; // check whether the well is operable under BHP limit with current reservoir condition virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) =0; diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 093a7304d..1391ef203 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -685,10 +685,9 @@ updateWellTestStateEconomic(const WellState& well_state, bool rate_limit_violated = false; const auto& quantity_limit = econ_production_limits.quantityLimit(); - const int np = number_of_phases_; if (econ_production_limits.onAnyRateLimit()) { if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN) - rate_limit_violated = checkRateEconLimits(econ_production_limits, &well_state.wellPotentials()[index_of_well_ * np], deferred_logger); + rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellPotentials(index_of_well_).data(), deferred_logger); else { rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger); } diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 18ce7033f..66f3596da 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -745,7 +745,7 @@ namespace Opm double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0); if (total_rate <= 0.0){ for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); double control_fraction = fractions[pu.phase_pos[Oil]]; if (control_fraction != 0.0) { for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); double control_fraction = fractions[pu.phase_pos[Water]]; if (control_fraction != 0.0) { for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); double control_fraction = fractions[pu.phase_pos[Gas]]; if (control_fraction != 0.0) { for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]]; if (control_fraction != 0.0) { for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); for (int p = 0; p fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials()); + const std::vector fractions = initialWellRateFractions(ebos_simulator, well_state); for (int p = 0; p std::vector WellInterface:: - initialWellRateFractions(const Simulator& ebosSimulator, const std::vector& potentials) const + initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const { const int np = this->number_of_phases_; std::vector scaling_factor(np); double total_potentials = 0.0; for (int p = 0; pindex_of_well_*np + p]; + total_potentials += well_state.wellPotentials(this->index_of_well_)[p]; } if (total_potentials > 0) { for (int p = 0; pindex_of_well_*np + p] / total_potentials; + scaling_factor[p] = well_state.wellPotentials(this->index_of_well_)[p] / total_potentials; } return scaling_factor; } diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index ff907f1a1..4277ec920 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -57,6 +57,9 @@ void WellState::base_init(const std::vector& cellPressures, this->thp_.clear(); this->temperature_.clear(); this->segment_state.clear(); + this->well_potentials_.clear(); + this->productivity_index_.clear(); + this->conn_productivity_index_.clear(); { // const int nw = wells->number_of_wells; const int nw = wells_ecl.size(); @@ -103,7 +106,7 @@ void WellState::initSingleWell(const std::vector& cellPressures, this->well_perf_data_.add(well.name(), well_perf_data); this->parallel_well_info_.add(well.name(), well_info); this->wellrates_.add(well.name(), std::vector(np, 0)); - + this->well_potentials_.add(well.name(), std::vector(np, 0)); const int num_perf_this_well = well_info->communication().sum(well_perf_data_[w].size()); this->segment_state.add(well.name(), SegmentState{}); this->perfpress_.add(well.name(), std::vector(num_perf_this_well, -1e100)); @@ -117,6 +120,8 @@ void WellState::initSingleWell(const std::vector& cellPressures, this->perfRateBrine_.add(well.name(), std::vector(num_perf_this_well, 0)); this->bhp_.add(well.name(), 0.0); this->thp_.add(well.name(), 0.0); + this->productivity_index_.add(well.name(), std::vector(np, 0)); + this->conn_productivity_index_.add(well.name(), std::vector(num_perf_this_well * np, 0)); if ( well.isInjector() ) this->temperature_.add(well.name(), well.injectionControls(summary_state).temperature); else @@ -341,9 +346,6 @@ void WellState::init(const std::vector& cellPressures, } } - productivity_index_.assign(nw * this->numPhases(), 0.0); - conn_productivity_index_.assign(nperf * this->numPhases(), 0.0); - well_potentials_.assign(nw * this->numPhases(), 0.0); for (int w = 0; w < nw; ++w) { switch (wells_ecl[w].getStatus()) { @@ -371,7 +373,6 @@ void WellState::init(const std::vector& cellPressures, continue; } const auto& wname = well.name(); - auto it = prevState->wellMap().find(well.name()); if (it != end) { @@ -403,9 +404,8 @@ void WellState::init(const std::vector& cellPressures, wellReservoirRates(w) = prevState->wellReservoirRates(oldIndex); // Well potentials - for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; iwellPotentials()[ oldidx ]; + for (int p=0; p < np; p++) { + this->wellPotentials(newIndex)[p] = prevState->wellPotentials(oldIndex)[p]; } // perfPhaseRates @@ -481,15 +481,8 @@ void WellState::init(const std::vector& cellPressures, } } - // Productivity index. - { - auto* thisWellPI = &this ->productivityIndex()[newIndex*np + 0]; - const auto* thatWellPI = &prevState->productivityIndex()[oldIndex*np + 0]; - - for (int p = 0; p < np; ++p) { - thisWellPI[p] = thatWellPI[p]; - } - } + // Productivity index. + this->productivity_index_.copy_welldata( prevState->productivity_index_, wname ); } // If in the new step, there is no THP related @@ -644,28 +637,26 @@ WellState::report(const int* globalCellIdxMap, } auto& well = res.at(wt.first); - const int well_rate_index = w * pu.num_phases; const auto& reservoir_rates = this->well_reservoir_rates_[w]; + const auto& well_potentials = this->well_potentials_[w]; + const auto& wpi = this->productivity_index_[w]; if (pu.phase_used[Water]) { - const auto i = well_rate_index + pu.phase_pos[Water]; well.rates.set(rt::reservoir_water, reservoir_rates[pu.phase_pos[Water]]); - well.rates.set(rt::productivity_index_water, this->productivity_index_[i]); - well.rates.set(rt::well_potential_water, this->well_potentials_[i]); + well.rates.set(rt::productivity_index_water, wpi[pu.phase_pos[Water]]); + well.rates.set(rt::well_potential_water, well_potentials[pu.phase_pos[Water]]); } if (pu.phase_used[Oil]) { - const auto i = well_rate_index + pu.phase_pos[Oil]; well.rates.set(rt::reservoir_oil, reservoir_rates[pu.phase_pos[Oil]]); - well.rates.set(rt::productivity_index_oil, this->productivity_index_[i]); - well.rates.set(rt::well_potential_oil, this->well_potentials_[i]); + well.rates.set(rt::productivity_index_oil, wpi[pu.phase_pos[Oil]]); + well.rates.set(rt::well_potential_oil, well_potentials[pu.phase_pos[Oil]]); } if (pu.phase_used[Gas]) { - const auto i = well_rate_index + pu.phase_pos[Gas]; well.rates.set(rt::reservoir_gas, reservoir_rates[pu.phase_pos[Gas]]); - well.rates.set(rt::productivity_index_gas, this->productivity_index_[i]); - well.rates.set(rt::well_potential_gas, this->well_potentials_[i]); + well.rates.set(rt::productivity_index_gas, wpi[pu.phase_pos[Gas]]); + well.rates.set(rt::well_potential_gas, well_potentials[pu.phase_pos[Gas]]); } if (pu.has_solvent || pu.has_zFraction) { @@ -752,14 +743,12 @@ void WellState::reportConnections(data::Well& well, 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(well_index)[np*local_comp_index]; - const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset; + const auto * rates = &this->perfPhaseRates(well_index)[np*local_comp_index]; + const auto& connPI = this->connectionProductivityIndex(well_index); for( int i = 0; i < np; ++i ) { comp.rates.set( phs[ i ], rates[i] ); - comp.rates.set( pi [ i ], *(connPI + i) ); + comp.rates.set( pi [ i ], connPI[i] ); } if ( pu.has_polymer ) { const auto& perf_polymer_rate = this->perfRatePolymer(well_index); @@ -794,7 +783,6 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, // what we do here, is to set the segment rates and perforation rates for (int w = 0; w < nw; ++w) { const auto& well_ecl = wells_ecl[w]; - const auto& wname = wells_ecl[w].name(); if ( well_ecl.isMultiSegment() ) { const WellSegments& segment_set = well_ecl.getSegments(); @@ -966,17 +954,15 @@ void WellState::shutWell(int well_index) this->wellrates_[well_index].assign(np, 0); auto& resv = this->well_reservoir_rates_[well_index]; - auto* wpi = &this->productivity_index_[np*well_index + 0]; + auto& wpi = this->productivity_index_[well_index]; for (int p = 0; p < np; ++p) { resv[p] = 0.0; wpi[p] = 0.0; } - const auto first = this->first_perf_index_[well_index]*np; - const auto last = first + this->numPerf(well_index)*np; - std::fill(this->conn_productivity_index_.begin() + first, - this->conn_productivity_index_.begin() + last, 0.0); + auto& connpi = this->conn_productivity_index_[well_index]; + connpi.assign(connpi.size(), 0); } void WellState::updateStatus(int well_index, Well::Status status) diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 541d47d5c..2856cecee 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -235,28 +235,29 @@ public: return this->segment_state[wname]; } - std::vector& productivityIndex() { - return productivity_index_; + std::vector& productivityIndex(std::size_t well_index) { + return this->productivity_index_[well_index]; } - const std::vector& productivityIndex() const { - return productivity_index_; + const std::vector& productivityIndex(std::size_t well_index) const { + return this->productivity_index_[well_index]; } - std::vector& connectionProductivityIndex() { - return this->conn_productivity_index_; + std::vector& connectionProductivityIndex(std::size_t well_index) { + return this->conn_productivity_index_[well_index]; } - const std::vector& connectionProductivityIndex() const { - return this->conn_productivity_index_; + const std::vector& connectionProductivityIndex(std::size_t well_index) const { + return this->conn_productivity_index_[well_index]; } - std::vector& wellPotentials() { - return well_potentials_; + + std::vector& wellPotentials(std::size_t well_index) { + return this->well_potentials_[well_index]; } - const std::vector& wellPotentials() const { - return well_potentials_; + const std::vector& wellPotentials(std::size_t well_index) const { + return this->well_potentials_[well_index]; } std::vector& perfThroughput(std::size_t well_index) { @@ -492,13 +493,13 @@ private: WellContainer segment_state; // Productivity Index - std::vector productivity_index_; + WellContainer> productivity_index_; // Connection-level Productivity Index - std::vector conn_productivity_index_; + WellContainer> conn_productivity_index_; // Well potentials - std::vector well_potentials_; + WellContainer> well_potentials_; data::Segment