diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh index f4e4e105b..e335c108a 100644 --- a/ebos/ecltracermodel.hh +++ b/ebos/ecltracermodel.hh @@ -433,7 +433,7 @@ protected: // TODO: Some inconsistencies here that perhaps should be clarified. The "offical" rate as reported below is // occasionally significant different from the sum over connections (as calculated above). Only observed // for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations. - std::size_t well_index = simulator_.problem().wellModel().wellState().wellIndex(well.name()); + std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value(); Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_]; rateWellTotal = official_well_rate_total; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 300122adc..9e575cd09 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -169,9 +169,9 @@ loadRestartData(const data::Wells& rst_wells, phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; } - for( const auto& wm : well_state.wellMap() ) { - const auto well_index = wm.second[ 0 ]; - const auto& rst_well = rst_wells.at( wm.first ); + for( std::size_t well_index = 0; well_index < well_state.size(); well_index++) { + const auto& well_name = well_state.name(well_index); + const auto& rst_well = rst_wells.at(well_name); auto& ws = well_state.well(well_index); ws.bhp = rst_well.bhp; ws.thp = rst_well.thp; @@ -206,7 +206,6 @@ loadRestartData(const data::Wells& rst_wells, if (handle_ms_well && !rst_well.segments.empty()) { // we need the well_ecl_ information - const std::string& well_name = wm.first; const Well& well_ecl = getWellEcl(well_name); const WellSegments& segment_set = well_ecl.getSegments(); @@ -1528,7 +1527,7 @@ updateAndCommunicateGroupData(const int reportStepIdx, // Set ALQ for off-process wells to zero for (const auto& wname : schedule().wellNames(reportStepIdx)) { const bool is_producer = schedule().getWell(wname, reportStepIdx).isProducer(); - const bool not_on_this_process = well_state.wellMap().count(wname) == 0; + const bool not_on_this_process = !well_state.has(wname); if (is_producer && not_on_this_process) { well_state.setALQ(wname, 0.0); } diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp index aadb2f0c2..06b9c54f0 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp @@ -1235,7 +1235,7 @@ bool GasLiftSingleWellGeneric::OptimizeState:: checkThpControl() { - const int well_index = this->parent.well_state_.wellIndex(this->parent.well_name_); + const int well_index = this->parent.well_state_.index(this->parent.well_name_).value(); const Well::ProducerCMode& control_mode = this->parent.well_state_.well(well_index).production_cmode; return control_mode == Well::ProducerCMode::THP; } diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index c73754929..305de3ff4 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -74,16 +74,13 @@ namespace { const auto& gefac = groupTmp.getGroupEfficiencyFactor(); rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector); } - const auto& end = wellState.wellMap().end(); for (const std::string& wellName : group.wells()) { - const auto& it = wellState.wellMap().find(wellName); - if (it == end) // the well is not found + const auto& well_index = wellState.index(wellName); + if (!well_index.has_value()) continue; - int well_index = it->second[0]; - - if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once + if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once { continue; } @@ -97,7 +94,7 @@ namespace { continue; double factor = wellEcl.getEfficiencyFactor(); - const auto& ws = wellState.well(well_index); + const auto& ws = wellState.well(well_index.value()); if (res_rates) { const auto& well_rates = ws.reservoir_rates; if (injector) @@ -216,15 +213,13 @@ namespace WellGroupHelpers const auto& gefac = groupTmp.getGroupEfficiencyFactor(); rate += gefac * sumSolventRates(groupTmp, schedule, wellState, reportStepIdx, injector); } - const auto& end = wellState.wellMap().end(); + for (const std::string& wellName : group.wells()) { - const auto& it = wellState.wellMap().find(wellName); - if (it == end) // the well is not found + const auto& well_index = wellState.index(wellName); + if (!well_index.has_value()) continue; - int well_index = it->second[0]; - - if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once + if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once { continue; } @@ -239,9 +234,9 @@ namespace WellGroupHelpers double factor = wellEcl.getEfficiencyFactor(); if (injector) - rate += factor * wellState.solventWellRate(well_index); + rate += factor * wellState.solventWellRate(well_index.value()); else - rate -= factor * wellState.solventWellRate(well_index); + rate -= factor * wellState.solventWellRate(well_index.value()); } return rate; } @@ -401,22 +396,19 @@ namespace WellGroupHelpers if (wellTmp.getStatus() == Well::Status::SHUT) continue; - const auto& end = wellState.wellMap().end(); - const auto& it = wellState.wellMap().find(wellName); - if (it == end) // the well is not found + const auto& well_index = wellState.index(wellName); + if (!well_index.has_value()) continue; - int well_index = it->second[0]; - - if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once + if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once { continue; } const double efficiency = wellTmp.getEfficiencyFactor(); // add contributino from wells not under group control - const auto& ws_nupcol = wellStateNupcol.well(well_index); - const auto& ws = wellState.well(well_index); + const auto& ws_nupcol = wellStateNupcol.well(well_index.value()); + const auto& ws = wellState.well(well_index.value()); if (isInjector) { if (ws.injection_cmode != Well::InjectorCMode::GRUP) for (int phase = 0; phase < np; phase++) { @@ -476,20 +468,17 @@ namespace WellGroupHelpers if (wellTmp.getStatus() == Well::Status::SHUT) continue; - const auto& end = wellState.wellMap().end(); - const auto& it = wellState.wellMap().find(wellName); - if (it == end) // the well is not found + const auto& well_index = wellState.index(wellName); + if (!well_index.has_value()) continue; - int well_index = it->second[0]; - - if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once + if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once { continue; } // scale rates - auto& ws = wellState.well(well_index); + auto& ws = wellState.well(well_index.value()); if (isInjector) { if (ws.injection_cmode == Well::InjectorCMode::GRUP) for (int phase = 0; phase < np; phase++) { @@ -568,19 +557,17 @@ namespace WellGroupHelpers updateWellRates(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState); } const int np = wellState.numPhases(); - const auto& end = wellState.wellMap().end(); for (const std::string& wellName : group.wells()) { std::vector rates(np, 0.0); - const auto& it = wellState.wellMap().find(wellName); - if (it != end) { // the well is found on this node - int well_index = it->second[0]; + const auto& well_index = wellState.index(wellName); + if (well_index.has_value()) { // the well is found on this node const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); int sign = 1; // production wellRates are negative. The users of currentWellRates uses the convention in // opm-common that production and injection rates are positive. if (!wellTmp.isInjector()) sign = -1; - const auto& ws = wellStateNupcol.well(well_index); + const auto& ws = wellStateNupcol.well(well_index.value()); for (int phase = 0; phase < np; ++phase) { rates[phase] = sign * ws.surface_rates[phase]; } @@ -1380,19 +1367,16 @@ namespace WellGroupHelpers if (wellTmp.getStatus() == Well::Status::SHUT) continue; - const auto& end = wellState.wellMap().end(); - const auto& it = wellState.wellMap().find(wellName); - if (it == end) // the well is not found + const auto& well_index = wellState.index(wellName); + if (!well_index.has_value()) // the well is not found continue; - int well_index = it->second[0]; - - if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once + if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once { continue; } - const auto& ws = wellState.well(well_index); + const auto& ws = wellState.well(well_index.value()); // add contribution from wells unconditionally for (int phase = 0; phase < np; phase++) { pot[phase] += wefac * ws.well_potentials[phase]; @@ -1430,18 +1414,16 @@ namespace WellGroupHelpers const Comm& comm, GuideRate* guideRate) { - const auto& end = wellState.wellMap().end(); for (const auto& well : schedule.getWells(reportStepIdx)) { double oilpot = 0.0; double gaspot = 0.0; double waterpot = 0.0; - const auto& it = wellState.wellMap().find(well.name()); - if (it != end && wellState.wellIsOwned(it->second[0], well.name())) + const auto& well_index = wellState.index(well.name()); + if (well_index.has_value() && wellState.wellIsOwned(well_index.value(), well.name())) { // the well is found and owned - int well_index = it->second[0]; - const auto& ws = wellState.well(well_index); + const auto& ws = wellState.well(well_index.value()); const auto& wpot = ws.well_potentials; if (pu.phase_used[BlackoilPhases::Liquid] > 0) oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]]; diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 9e8819fb4..5e7d7a9f5 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -39,29 +39,17 @@ void WellState::base_init(const std::vector& cellPressures, const SummaryState& summary_state) { // clear old name mapping - this->wellMap_.clear(); this->parallel_well_info_.clear(); this->wells_.clear(); { // const int nw = wells->number_of_wells; const int nw = wells_ecl.size(); // const int np = wells->number_of_phases; - int connpos = 0; for (int w = 0; w < nw; ++w) { const Well& well = wells_ecl[w]; // Initialize bhp(), thp(), wellRates(), temperature(). initSingleWell(cellPressures, well, well_perf_data[w], parallel_well_info[w], summary_state); - - // Setup wellname -> well index mapping. - const int num_perf_this_well = well_perf_data[w].size(); - std::string name = well.name(); - assert( !name.empty() ); - mapentry_t& wellMapEntry = wellMap_[name]; - wellMapEntry[ 0 ] = w; - wellMapEntry[ 1 ] = connpos; - wellMapEntry[ 2 ] = num_perf_this_well; - connpos += num_perf_this_well; } } } @@ -85,10 +73,9 @@ void WellState::initSingleWell(const std::vector& cellPressures, double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56; this->parallel_well_info_.add(well.name(), well_info); - const int num_perf_this_well = well_info->communication().sum(well_perf_data.size()); - auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), static_cast(num_perf_this_well), static_cast(np), temp}); + auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), well_perf_data.size(), static_cast(np), temp}); - if ( num_perf_this_well == 0 ) + if ( ws.perf_data.empty()) return; const auto inj_controls = well.isInjector() ? well.injectionControls(summary_state) : Well::InjectionControls(0); @@ -256,11 +243,10 @@ void WellState::init(const std::vector& cellPressures, // rates divided by the number of perforations. const auto& ecl_well = wells_ecl[w]; const auto& wname = ecl_well.name(); - const auto& well_info = this->wellMap().at(wname); - const int num_perf_this_well = well_info[2]; - const int global_num_perf_this_well = ecl_well.getConnections().num_open(); auto& ws = this->well(w); auto& perf_data = ws.perf_data; + const int num_perf_this_well = perf_data.size(); + const int global_num_perf_this_well = ecl_well.getConnections().num_open(); const auto& perf_input = well_perf_data[w]; for (int perf = 0; perf < num_perf_this_well; ++perf) { @@ -309,20 +295,16 @@ void WellState::init(const std::vector& cellPressures, // intialize wells that have been there before // order may change so the mapping is based on the well name - if (prevState && !prevState->wellMap().empty()) { - auto end = prevState->wellMap().end(); + if (prevState && prevState->size() > 0) { for (int w = 0; w < nw; ++w) { const Well& well = wells_ecl[w]; if (well.getStatus() == Well::Status::SHUT) { continue; } auto& new_well = this->well(w); - auto it = prevState->wellMap().find(well.name()); - if (it != end) - { - const int oldIndex = it->second[ 0 ]; - - const auto& prev_well = prevState->well(oldIndex); + const auto& old_index = prevState->index(well.name()); + if (old_index.has_value()) { + const auto& prev_well = prevState->well(old_index.value()); new_well.init_timestep(prev_well); @@ -346,16 +328,8 @@ void WellState::init(const std::vector& cellPressures, new_well.well_potentials = prev_well.well_potentials; // perfPhaseRates - const int num_perf_old_well = (*it).second[ 2 ]; - const auto new_iter = this->wellMap().find(well.name()); - if (new_iter == this->wellMap().end()) { - throw std::logic_error { - well.name() + " is not in internal well map - " - "Bug in WellState" - }; - } - - const int num_perf_this_well = new_iter->second[2]; + const int num_perf_old_well = prev_well.perf_data.size(); + const int num_perf_this_well = new_well.perf_data.size(); const bool global_num_perf_same = (num_perf_this_well == num_perf_old_well); // copy perforation rates when the number of @@ -459,8 +433,7 @@ WellState::report(const int* globalCellIdxMap, const auto& pu = this->phaseUsage(); data::Wells res; - for( const auto& [wname, winfo]: this->wellMap() ) { - const auto well_index = winfo[ 0 ]; + for( std::size_t well_index = 0; well_index < this->size(); well_index++) { const auto& ws = this->well(well_index); if ((ws.status == Well::Status::SHUT) && !wasDynamicallyClosed(well_index)) { @@ -471,6 +444,8 @@ WellState::report(const int* globalCellIdxMap, const auto& well_potentials = ws.well_potentials; const auto& wpi = ws.productivity_index; const auto& wv = ws.surface_rates; + const auto& wname = this->name(well_index); + data::Well well; well.bhp = ws.bhp; @@ -572,8 +547,6 @@ void WellState::reportConnections(std::vector& connections, connection.reservoir_rate = perf_rates[i]; connection.trans_factor = perf_data.connection_transmissibility_factor[i]; } - assert(num_perf_well == int(connections.size())); - const int np = pu.num_phases; size_t local_comp_index = 0; @@ -616,7 +589,6 @@ void WellState::reportConnections(std::vector& connections, ++local_comp_index; } - assert(local_comp_index == this->perfdata[well_index].size()); } void WellState::initWellStateMSWell(const std::vector& wells_ecl, @@ -672,28 +644,22 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, auto& perf_data = ws.perf_data; // for the seg_rates_, now it becomes a recursive solution procedure. - { - // make sure the information from wells_ecl consistent with wells - assert((n_activeperf == this->wellMap().at(well_ecl.name())[2]) && - "Inconsistent number of reservoir connections in well"); - - if (pu.phase_used[Gas]) { - auto& perf_rates = perf_data.phase_rates; - const int gaspos = pu.phase_pos[Gas]; - // scale the phase rates for Gas to avoid too bad initial guess for gas fraction - // it will probably benefit the standard well too, while it needs to be justified - // TODO: to see if this strategy can benefit StandardWell too - // TODO: it might cause big problem for gas rate control or if there is a gas rate limit - // maybe the best way is to initialize the fractions first then get the rates - for (int perf = 0; perf < n_activeperf; perf++) - perf_rates[perf*np + gaspos] *= 100; - } - - const auto& perf_rates = perf_data.phase_rates; - std::vector perforation_rates(perf_rates.begin(), perf_rates.end()); - - calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates); + if (pu.phase_used[Gas]) { + auto& perf_rates = perf_data.phase_rates; + const int gaspos = pu.phase_pos[Gas]; + // scale the phase rates for Gas to avoid too bad initial guess for gas fraction + // it will probably benefit the standard well too, while it needs to be justified + // TODO: to see if this strategy can benefit StandardWell too + // TODO: it might cause big problem for gas rate control or if there is a gas rate limit + // maybe the best way is to initialize the fractions first then get the rates + for (int perf = 0; perf < n_activeperf; perf++) + perf_rates[perf*np + gaspos] *= 100; } + + const auto& perf_rates = perf_data.phase_rates; + std::vector perforation_rates(perf_rates.begin(), perf_rates.end()); + + calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates); // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment // if there is no perforation associated with this segment, it uses the pressure from the outlet segment // which requres the ordering is successful @@ -947,21 +913,11 @@ bool WellState::wellIsOwned(std::size_t well_index, bool WellState::wellIsOwned(const std::string& wellName) const { - const auto& it = this->wellMap_.find( wellName ); - if (it == this->wellMap_.end()) { + const auto& well_index = this->index(wellName); + if (!well_index.has_value()) OPM_THROW(std::logic_error, "Could not find well " << wellName << " in well map"); - } - const int well_index = it->second[0]; - return wellIsOwned(well_index, wellName); -} -int WellState::wellIndex(const std::string& wellName) const -{ - const auto& it = this->wellMap_.find( wellName ); - if (it == this->wellMap_.end()) { - OPM_THROW(std::logic_error, "Could not find well " << wellName << " in well map"); - } - return it->second[0]; + return wellIsOwned(well_index.value(), wellName); } int WellState::numSegments(const int well_id) const diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index f05f8dcf9..73a5c131a 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -52,9 +52,6 @@ class Schedule; class WellState { public: - using mapentry_t = std::array; - using WellMapType = std::map; - static const uint64_t event_mask = ScheduleEvents::WELL_STATUS_CHANGE + ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::INJECTION_UPDATE; virtual ~WellState() = default; @@ -69,11 +66,8 @@ public: this->phase_usage_ = pu; } - const WellMapType& wellMap() const { return wellMap_; } - WellMapType& wellMap() { return wellMap_; } - std::size_t size() const { - return this->wellMap_.size(); + return this->wells_.size(); } @@ -82,8 +76,6 @@ public: return this->size(); } - int wellIndex(const std::string& wellName) const; - const ParallelWellInfo& parallelWellInfo(std::size_t well_index) const; @@ -250,6 +242,10 @@ public: return this->wells_.well_name(well_index); } + std::optional index(const std::string& well_name) const { + return this->wells_.well_index(well_name); + } + const SingleWellState& well(std::size_t well_index) const { return this->wells_[well_index]; } @@ -271,7 +267,6 @@ public: } private: - WellMapType wellMap_; // Use of std::optional<> here is a technical crutch, the // WellStateFullyImplicitBlackoil class should be default constructible, // whereas the GlobalWellInfo is not.