From 1a6737aff399a6e0c795d3bbfa26f869282b6fcb Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Thu, 3 Jun 2021 11:08:55 +0200 Subject: [PATCH] Use well index for perforation rates --- .../wells/BlackoilWellModel_impl.hpp | 4 ++-- .../wells/MultisegmentWell_impl.hpp | 3 +-- opm/simulators/wells/StandardWell_impl.hpp | 5 ++--- .../wells/WellInterfaceFluidSystem.cpp | 2 +- opm/simulators/wells/WellState.cpp | 22 ++++++++----------- opm/simulators/wells/WellState.hpp | 9 ++++++-- 6 files changed, 22 insertions(+), 23 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 80985ec7d..c219c1e73 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1997,7 +1997,7 @@ namespace Opm { auto& perf_pressure = well_state.perfPress(well_index); auto& perf_rates = well_state.perfRates(well_index); - auto * perf_phase_rates = &well_state.mutable_perfPhaseRates()[wm.second[1]*np]; + auto * perf_phase_rates = well_state.perfPhaseRates(well_index); const auto& perf_data = this->well_perf_data_[well_index]; for (std::size_t perf_index = 0; perf_index < perf_data.size(); perf_index++) { @@ -3374,7 +3374,7 @@ namespace Opm { auto& well_info = *local_parallel_well_info_[wellID]; const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size()); - auto * perf_phase_rate = &this->wellState().perfPhaseRates()[connpos]; + auto * perf_phase_rate = this->wellState().perfPhaseRates(wellID); for (int perf = 0; perf < num_perf_this_well; ++perf) { const int cell_idx = well_perf_data_[wellID][perf].cell_index; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 9dca12e19..569815636 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -2637,8 +2637,7 @@ namespace Opm // calculating the perforation rate for each perforation that belongs to this segment const EvalWell seg_pressure = getSegmentPressure(seg); - const int rate_start_offset = first_perf_ * number_of_phases_; - auto * perf_rates = &well_state.mutable_perfPhaseRates()[rate_start_offset]; + auto * perf_rates = well_state.perfPhaseRates(this->index_of_well_); auto& perf_press_state = well_state.perfPress(this->index_of_well_); for (const int perf : segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 436667419..c8a22ab89 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -591,8 +591,7 @@ namespace Opm const int np = number_of_phases_; std::vector connectionRates = connectionRates_; // Copy to get right size. - const int rate_start_offset = first_perf_ * number_of_phases_; - auto * perf_rates = &well_state.mutable_perfPhaseRates()[rate_start_offset]; + auto * perf_rates = well_state.perfPhaseRates(this->index_of_well_); for (int perf = 0; perf < number_of_perforations_; ++perf) { // Calculate perforation quantities. std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.0}); @@ -2105,7 +2104,7 @@ namespace Opm const int nperf = number_of_perforations_; const int np = number_of_phases_; std::vector perfRates(b_perf.size(),0.0); - const auto * perf_rates_state = &well_state.perfPhaseRates()[first_perf_ * np]; + const auto * perf_rates_state = well_state.perfPhaseRates(this->index_of_well_); for (int perf = 0; perf < nperf; ++perf) { for (int comp = 0; comp < np; ++comp) { diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index fd8e63614..7b735fdb2 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -846,7 +846,7 @@ checkMaxRatioLimitCompletions(const WellState& well_state, double max_ratio_completion = 0; const int np = number_of_phases_; - const auto * perf_phase_rates = &well_state.perfPhaseRates()[first_perf_ * np]; + const auto * perf_phase_rates = well_state.perfPhaseRates(this->index_of_well_); // look for the worst_offending_completion for (const auto& completion : completions_) { std::vector completion_rates(np, 0.0); diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 2e2e5fc35..8cd13d0de 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -298,7 +298,9 @@ void WellState::init(const std::vector& cellPressures, const int num_perf_this_well = well_info[2]; const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well); auto& perf_press = this->perfPress(w); - auto * phase_rates = &this->mutable_perfPhaseRates()[connpos * this->numPhases()]; + + first_perf_index_[w] = connpos; + auto phase_rates = this->perfPhaseRates(w); for (int perf = 0; perf < num_perf_this_well; ++perf) { if (wells_ecl[w].getStatus() == Well::Status::OPEN) { @@ -308,7 +310,6 @@ void WellState::init(const std::vector& cellPressures, } perf_press[perf] = cellPressures[well_perf_data[w][perf].cell_index]; } - first_perf_index_[w] = connpos; this->well_reservoir_rates_.add(wname, std::vector(np, 0)); this->well_dissolved_gas_rates_.add(wname, 0); @@ -408,7 +409,6 @@ void WellState::init(const std::vector& cellPressures, } // perfPhaseRates - const int oldPerf_idx_beg = (*it).second[ 1 ]; const int num_perf_old_well = (*it).second[ 2 ]; const auto new_iter = this->wellMap().find(well.name()); if (new_iter == this->wellMap().end()) { @@ -418,7 +418,6 @@ void WellState::init(const std::vector& cellPressures, }; } - const int connpos = new_iter->second[1]; const int num_perf_this_well = new_iter->second[2]; const int num_perf_changed = parallel_well_info[w]->communication() @@ -432,8 +431,8 @@ void WellState::init(const std::vector& cellPressures, // number of perforations. if (global_num_perf_same) { - const auto * src_rates = &prevState->perfPhaseRates()[oldPerf_idx_beg* np]; - auto * target_rates = &this->mutable_perfPhaseRates()[connpos*np]; + const auto * src_rates = prevState->perfPhaseRates(oldIndex); + auto * target_rates = this->perfPhaseRates(newIndex); for (int perf_index = 0; perf_index < num_perf_this_well; perf_index++) { for (int p = 0; p < np; p++) { target_rates[perf_index*np + p] = src_rates[perf_index*np + p]; @@ -441,7 +440,7 @@ void WellState::init(const std::vector& cellPressures, } } else { const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well); - auto * target_rates = &this->mutable_perfPhaseRates()[connpos*np]; + auto * target_rates = this->perfPhaseRates(newIndex); for (int perf_index = 0; perf_index < num_perf_this_well; perf_index++) { for (int p = 0; p < np; ++p) { target_rates[perf_index*np + p] = wellRates(w)[p] / double(global_num_perf_this_well); @@ -771,7 +770,7 @@ void WellState::reportConnections(data::Well& well, for( auto& comp : well.connections) { const auto connPhaseOffset = np * (wt.second[1] + local_comp_index); - const auto * rates = &this->perfPhaseRates()[connPhaseOffset]; + const auto * rates = &this->perfPhaseRates(well_index)[np*local_comp_index]; const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset; for( int i = 0; i < np; ++i ) { @@ -813,7 +812,6 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, const auto& well_ecl = wells_ecl[w]; const auto& wname = wells_ecl[w].name(); const auto& well_info = this->wellMap().at(wname); - const int connpos = well_info[1]; const int num_perf_this_well = well_info[2]; if ( well_ecl.isMultiSegment() ) { @@ -852,14 +850,12 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, // for the seg_rates_, now it becomes a recursive solution procedure. { - const int start_perf = connpos; - // make sure the information from wells_ecl consistent with wells assert((n_activeperf == num_perf_this_well) && "Inconsistent number of reservoir connections in well"); if (pu.phase_used[Gas]) { - auto * perf_rates = &this->mutable_perfPhaseRates()[np * start_perf]; + auto * perf_rates = this->perfPhaseRates(w); 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 @@ -870,7 +866,7 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, perf_rates[perf*np + gaspos] *= 100; } - const auto * perf_rates = &perfPhaseRates()[np*start_perf]; + const auto * perf_rates = this->perfPhaseRates(w); std::vector perforation_rates(perf_rates, perf_rates + num_perf_this_well*np); auto& segments = this->segments(w); diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index a071992ed..d26aecb62 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -100,8 +100,13 @@ public: const SummaryState& summary_state); /// One rate per phase and well connection. - std::vector& mutable_perfPhaseRates() { return perfphaserates_; } - const std::vector& perfPhaseRates() const { return perfphaserates_; } + double * perfPhaseRates(std::size_t well_index) { + return &this->perfphaserates_[this->first_perf_index_[well_index] * this->numPhases()]; + } + + const double * perfPhaseRates(std::size_t well_index) const { + return &this->perfphaserates_[this->first_perf_index_[well_index] * this->numPhases()]; + } /// One current control per injecting well. Well::InjectorCMode currentInjectionControl(std::size_t well_index) const { return current_injection_controls_[well_index]; }