diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 7bd59cfb2..89b050596 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -2642,6 +2642,9 @@ 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_press_state = &well_state.perfPress()[first_perf_]; for (const int perf : segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); @@ -2662,11 +2665,10 @@ namespace Opm } // store the perf pressure and rates - const int rate_start_offset = (first_perf_ + perf) * number_of_phases_; for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - well_state.mutable_perfPhaseRates()[rate_start_offset + ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value(); + perf_rates[perf*number_of_phases_ + ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value(); } - well_state.perfPress()[first_perf_ + perf] = perf_press.value(); + perf_press_state[perf] = perf_press.value(); for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { // the cq_s entering mass balance equations need to consider the efficiency factors. diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 36147476d..350eb8d40 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -584,6 +584,8 @@ 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]; for (int perf = 0; perf < number_of_perforations_; ++perf) { // Calculate perforation quantities. std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.0}); @@ -620,7 +622,7 @@ namespace Opm if (has_solvent && componentIdx == contiSolventEqIdx) { well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); } else { - well_state.mutable_perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value(); + perf_rates[perf*np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value(); } } @@ -836,7 +838,8 @@ namespace Opm } // Store the perforation pressure for later usage. - well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; + auto * perf_press = &well_state.perfPress()[first_perf_]; + perf_press[perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; } @@ -1644,9 +1647,10 @@ namespace Opm } // Compute the average pressure in each well block + const auto * perf_press = &well_state.perfPress()[first_perf_]; auto p_above = this->parallel_well_info_.communicateAboveValues(well_state.bhp()[w], - well_state.perfPress().data()+first_perf_, - nperf); + perf_press, + nperf); for (int perf = 0; perf < nperf; ++perf) { const int cell_idx = well_cells_[perf]; @@ -1655,7 +1659,7 @@ namespace Opm // 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 = (well_state.perfPress()[first_perf_ + perf] + p_above[perf])/2; + 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(); @@ -2080,9 +2084,10 @@ 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]; for (int perf = 0; perf < nperf; ++perf) { for (int comp = 0; comp < np; ++comp) { - perfRates[perf * num_components_ + comp] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(comp)]; + perfRates[perf * num_components_ + comp] = perf_rates_state[perf * np + ebosCompIdxToFlowCompIdx(comp)]; } if(has_solvent) { perfRates[perf * num_components_ + contiSolventEqIdx] = well_state.perfRateSolvent()[first_perf_ + perf]; diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 5acb4b730..4b8846227 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -322,12 +322,14 @@ namespace Opm const int num_perf_well = pd.size(); well.connections.resize(num_perf_well); + const auto * perf_rates = &this->perfRates()[itr.second[1]]; + const auto * perf_pressure = &this->perfPress()[itr.second[1]]; 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.pressure = perf_pressure[i]; + connection.reservoir_rate = perf_rates[i]; connection.trans_factor = pd[i].connection_transmissibility_factor; } assert(num_perf_well == int(well.connections.size())); diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 5dcd5667a..206867e73 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -150,14 +150,16 @@ namespace Opm const int connpos = well_info[1]; 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()[connpos]; + auto * phase_rates = &this->mutable_perfPhaseRates()[connpos * this->numPhases()]; - for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) { + for (int perf = 0; perf < num_perf_this_well; ++perf) { if (wells_ecl[w].getStatus() == Well::Status::OPEN) { for (int p = 0; p < this->numPhases(); ++p) { - perfphaserates_[this->numPhases()*perf + p] = wellRates()[this->numPhases()*w + p] / double(global_num_perf_this_well); + phase_rates[this->numPhases()*perf + p] = wellRates()[this->numPhases()*w + p] / double(global_num_perf_this_well); } } - perfPress()[perf] = cellPressures[well_perf_data[w][perf-connpos].cell_index]; + perf_press[perf] = cellPressures[well_perf_data[w][perf].cell_index]; } num_perf_[w] = num_perf_this_well; first_perf_index_[w] = connpos; @@ -289,17 +291,19 @@ namespace Opm // number of perforations. if (global_num_perf_same) { - int old_perf_phase_idx = oldPerf_idx_beg *np; - for (int perf_phase_idx = connpos*np; - perf_phase_idx < (connpos + num_perf_this_well)*np; ++perf_phase_idx, ++old_perf_phase_idx ) - { - mutable_perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ]; + const auto * src_rates = &prevState->perfPhaseRates()[oldPerf_idx_beg* np]; + auto * target_rates = &this->mutable_perfPhaseRates()[connpos*np]; + 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]; + } } } else { const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well); - for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) { + auto * target_rates = &this->mutable_perfPhaseRates()[connpos*np]; + for (int perf_index = 0; perf_index < num_perf_this_well; perf_index++) { for (int p = 0; p < np; ++p) { - mutable_perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(global_num_perf_this_well); + target_rates[perf_index*np + p] = wellRates()[np*newIndex + p] / double(global_num_perf_this_well); } } } @@ -307,10 +311,11 @@ namespace Opm // perfPressures if (global_num_perf_same) { - int oldPerf_idx = oldPerf_idx_beg; - for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx ) + auto * target_press = &perfPress()[connpos]; + const auto * src_press = &prevState->perfPress()[oldPerf_idx_beg]; + for (int perf = 0; perf < num_perf_this_well; ++perf) { - perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ]; + target_press[perf] = src_press[perf]; } } @@ -676,27 +681,26 @@ namespace Opm // for the seg_rates_, now it becomes a recursive solution procedure. { const int start_perf = connpos; - const int start_perf_next_well = connpos + num_perf_this_well; // make sure the information from wells_ecl consistent with wells assert(n_activeperf == (start_perf_next_well - start_perf)); if (pu.phase_used[Gas]) { + auto * perf_rates = &this->mutable_perfPhaseRates()[np * start_perf]; 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++) { - const int perf_pos = start_perf + perf; - mutable_perfPhaseRates()[np * perf_pos + gaspos] *= 100.; - } + for (int perf = 0; perf < n_activeperf; perf++) + perf_rates[perf*np + gaspos] *= 100; } - const std::vector perforation_rates(perfPhaseRates().begin() + np * start_perf, - perfPhaseRates().begin() + np * start_perf_next_well); // the perforation rates for this well + const auto * perf_rates = &perfPhaseRates()[np*start_perf]; + std::vector perforation_rates(perf_rates, perf_rates + num_perf_this_well*np); std::vector segment_rates; + calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, segment_rates); std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(seg_rates_)); } @@ -711,10 +715,11 @@ namespace Opm seg_press_.push_back(bhp()[w]); const int top_segment = top_segment_index_[w]; const int start_perf = connpos; + const auto * perf_press = &this->perfPress()[start_perf]; for (int seg = 1; seg < well_nseg; ++seg) { if ( !segment_perforations[seg].empty() ) { const int first_perf = segment_perforations[seg][0]; - seg_press_.push_back(perfPress()[start_perf + first_perf]); + seg_press_.push_back(perf_press[first_perf]); } else { // seg_press_.push_back(bhp); // may not be a good decision // using the outlet segment pressure // it needs the ordering is correct