Merge pull request #3339 from joakim-hove/perf-well-access

Perf well access
This commit is contained in:
Joakim Hove
2021-06-04 08:41:47 +02:00
committed by GitHub
6 changed files with 83 additions and 68 deletions

View File

@@ -1996,7 +1996,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++) {
@@ -3373,7 +3373,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;

View File

@@ -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];

View File

@@ -591,8 +591,7 @@ namespace Opm
const int np = number_of_phases_;
std::vector<RateVector> 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<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.0});
@@ -629,7 +628,7 @@ namespace Opm
// Store the perforation phase flux for later usage.
if (has_solvent && componentIdx == contiSolventEqIdx) {
auto * perf_rate_solvent = &well_state.perfRateSolvent()[first_perf_];
auto * perf_rate_solvent = well_state.perfRateSolvent(this->index_of_well_);
perf_rate_solvent[perf] = cq_s[componentIdx].value();
} else {
perf_rates[perf*np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
@@ -795,7 +794,7 @@ namespace Opm
cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
}
// Note. Efficiency factor is handled in the output layer
auto * perf_rate_polymer = &well_state.perfRatePolymer()[first_perf_];
auto * perf_rate_polymer = well_state.perfRatePolymer(this->index_of_well_);
perf_rate_polymer[perf] = cq_s_poly.value();
cq_s_poly *= well_efficiency_factor_;
@@ -828,7 +827,7 @@ namespace Opm
const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
cq_s_zfrac_effective *= extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
}
auto * perf_rate_solvent = &well_state.perfRateSolvent()[first_perf_];
auto * perf_rate_solvent = well_state.perfRateSolvent(this->index_of_well_);
perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
cq_s_zfrac_effective *= well_efficiency_factor_;
@@ -845,7 +844,7 @@ namespace Opm
cq_s_sm *= extendEval(intQuants.fluidState().saltConcentration());
}
// Note. Efficiency factor is handled in the output layer
auto * perf_rate_brine = &well_state.perfRateBrine()[this->first_perf_];
auto * perf_rate_brine = well_state.perfRateBrine(this->index_of_well_);
perf_rate_brine[perf] = cq_s_sm.value();
cq_s_sm *= well_efficiency_factor_;
@@ -1305,8 +1304,8 @@ namespace Opm
// other primary variables related to polymer injectivity study
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
auto * perf_water_velocity = &well_state.perfWaterVelocity()[this->first_perf_];
auto * perf_skin_pressure = &well_state.perfSkinPressure()[this->first_perf_];
auto * perf_water_velocity = well_state.perfWaterVelocity(this->index_of_well_);
auto * perf_skin_pressure = well_state.perfSkinPressure(this->index_of_well_);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
perf_water_velocity[perf] = primary_variables_[Bhp + 1 + perf];
perf_skin_pressure[perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
@@ -2105,7 +2104,7 @@ namespace Opm
const int nperf = number_of_perforations_;
const int np = number_of_phases_;
std::vector<double> 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) {
@@ -2114,7 +2113,7 @@ namespace Opm
}
if constexpr (has_solvent) {
const auto * solvent_perf_rates_state = &well_state.perfRateSolvent()[this->first_perf_];
const auto * solvent_perf_rates_state = well_state.perfRateSolvent(this->index_of_well_);
for (int perf = 0; perf < nperf; ++perf) {
perfRates[perf * num_components_ + contiSolventEqIdx] = solvent_perf_rates_state[perf];
}
@@ -2833,8 +2832,8 @@ namespace Opm
// other primary variables related to polymer injection
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
const auto * water_velocity = &well_state.perfWaterVelocity()[first_perf_];
const auto * skin_pressure = &well_state.perfSkinPressure()[first_perf_];
const auto * water_velocity = well_state.perfWaterVelocity(this->index_of_well_);
const auto * skin_pressure = well_state.perfSkinPressure(this->index_of_well_);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = skin_pressure[perf];
@@ -3202,7 +3201,7 @@ namespace Opm
{
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
auto * perf_water_throughput = &well_state.perfThroughput()[first_perf_];
auto * perf_water_throughput = well_state.perfThroughput(this->index_of_well_);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
@@ -3263,7 +3262,7 @@ namespace Opm
const EvalWell eq_wat_vel = primary_variables_evaluation_[wat_vel_index] - water_velocity;
resWell_[0][wat_vel_index] = eq_wat_vel.value();
const auto * perf_water_throughput = &well_state.perfThroughput()[this->first_perf_];
const auto * perf_water_throughput = well_state.perfThroughput(this->index_of_well_);
const double throughput = perf_water_throughput[perf];
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
@@ -3442,7 +3441,7 @@ namespace Opm
const int wat_vel_index = Bhp + 1 + perf;
const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index];
if (water_velocity > 0.) { // injecting
const auto * perf_water_throughput = &well_state.perfThroughput()[this->first_perf_];
const auto * perf_water_throughput = well_state.perfThroughput(this->index_of_well_);
const double throughput = perf_water_throughput[perf];
const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
cq_s_polymw *= molecular_weight;

View File

@@ -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<double> completion_rates(np, 0.0);

View File

@@ -298,7 +298,9 @@ void WellState::init(const std::vector<double>& 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<double>& 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<double>(np, 0));
this->well_dissolved_gas_rates_.add(wname, 0);
@@ -408,7 +409,6 @@ void WellState::init(const std::vector<double>& 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<double>& 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<double>& 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<double>& 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);
@@ -464,10 +463,11 @@ void WellState::init(const std::vector<double>& cellPressures,
if (pu.has_solvent) {
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 * solvent_target = this->perfRateSolvent(newIndex);
const auto * solvent_src = prevState->perfRateSolvent(oldIndex);
for (int perf = 0; perf < num_perf_this_well; ++perf)
{
perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ];
solvent_target[perf] = solvent_src[perf];
}
}
}
@@ -480,13 +480,13 @@ void WellState::init(const std::vector<double>& cellPressures,
if (pu.has_polymermw) {
if (global_num_perf_same)
{
auto * throughput_target = &perf_water_throughput_[connpos];
auto * pressure_target = &perf_skin_pressure_[connpos];
auto * velocity_target = &perf_water_velocity_[connpos];
auto * throughput_target = this->perfThroughput(newIndex);
auto * pressure_target = this->perfSkinPressure(newIndex);
auto * velocity_target = this->perfWaterVelocity(newIndex);
const auto * throughput_src = &prevState->perfThroughput()[oldPerf_idx_beg];
const auto * pressure_src = &prevState->perfSkinPressure()[oldPerf_idx_beg];
const auto * velocity_src = &prevState->perfWaterVelocity()[oldPerf_idx_beg];
const auto * throughput_src = prevState->perfThroughput(oldIndex);
const auto * pressure_src = prevState->perfSkinPressure(oldIndex);
const auto * velocity_src = prevState->perfWaterVelocity(oldIndex);
for (int perf = 0; perf < num_perf_this_well; ++perf)
{
@@ -770,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 ) {
@@ -778,15 +778,15 @@ void WellState::reportConnections(data::Well& well,
comp.rates.set( pi [ i ], *(connPI + i) );
}
if ( pu.has_polymer ) {
const auto * perf_polymer_rate = &this->perfRatePolymer()[wt.second[1]];
const auto * perf_polymer_rate = this->perfRatePolymer(well_index);
comp.rates.set( rt::polymer, perf_polymer_rate[local_comp_index]);
}
if ( pu.has_brine ) {
const auto * perf_brine_rate = &this->perfRateBrine()[wt.second[1]];
const auto * perf_brine_rate = this->perfRateBrine(well_index);
comp.rates.set( rt::brine, perf_brine_rate[local_comp_index]);
}
if ( pu.has_solvent ) {
const auto * perf_solvent_rate = &this->perfRateSolvent()[wt.second[1]];
const auto * perf_solvent_rate = this->perfRateSolvent(well_index);
comp.rates.set( rt::solvent, perf_solvent_rate[local_comp_index] );
}
@@ -812,7 +812,6 @@ void WellState::initWellStateMSWell(const std::vector<Well>& 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() ) {
@@ -851,14 +850,12 @@ void WellState::initWellStateMSWell(const std::vector<Well>& 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
@@ -869,7 +866,7 @@ void WellState::initWellStateMSWell(const std::vector<Well>& 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<double> perforation_rates(perf_rates, perf_rates + num_perf_this_well*np);
auto& segments = this->segments(w);
@@ -955,19 +952,19 @@ WellState::calculateSegmentRates(const std::vector<std::vector<int>>& segment_in
double WellState::solventWellRate(const int w) const
{
const auto * perf_rates_solvent = &perfRateSolvent_[first_perf_index_[w]];
const auto * perf_rates_solvent = this->perfRateSolvent(w);
return parallel_well_info_[w]->sumPerfValues(perf_rates_solvent, perf_rates_solvent + this->numPerf(w));
}
double WellState::polymerWellRate(const int w) const
{
const auto * perf_rates_polymer = &perfRatePolymer_[first_perf_index_[w]];
const auto * perf_rates_polymer = this->perfRatePolymer(w);
return parallel_well_info_[w]->sumPerfValues(perf_rates_polymer, perf_rates_polymer + this->numPerf(w));
}
double WellState::brineWellRate(const int w) const
{
const auto * perf_rates_brine = &perfRateBrine_[first_perf_index_[w]];
const auto * perf_rates_brine = this->perfRateBrine(w);
return parallel_well_info_[w]->sumPerfValues(perf_rates_brine, perf_rates_brine + this->numPerf(w));
}

View File

@@ -100,8 +100,13 @@ public:
const SummaryState& summary_state);
/// One rate per phase and well connection.
std::vector<double>& mutable_perfPhaseRates() { return perfphaserates_; }
const std::vector<double>& 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]; }
@@ -153,22 +158,37 @@ public:
}
/// One rate pr well connection.
std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
double * perfRateSolvent(std::size_t well_index) {
return &perfRateSolvent_[this->first_perf_index_[well_index]];
}
const double * perfRateSolvent(std::size_t well_index) const {
return &perfRateSolvent_[this->first_perf_index_[well_index]];
}
/// One rate pr well
double solventWellRate(const int w) const;
/// One rate pr well connection.
std::vector<double>& perfRatePolymer() { return perfRatePolymer_; }
const std::vector<double>& perfRatePolymer() const { return perfRatePolymer_; }
double * perfRatePolymer(std::size_t well_index) {
return &this->perfRatePolymer_[this->first_perf_index_[well_index]];
}
const double * perfRatePolymer(std::size_t well_index) const {
return &this->perfRatePolymer_[this->first_perf_index_[well_index]];
}
/// One rate pr well
double polymerWellRate(const int w) const;
/// One rate pr well connection.
std::vector<double>& perfRateBrine() { return perfRateBrine_; }
const std::vector<double>& perfRateBrine() const { return perfRateBrine_; }
double* perfRateBrine(std::size_t well_index) {
return &this->perfRateBrine_[this->first_perf_index_[well_index]];
}
const double* perfRateBrine(std::size_t well_index) const {
return &this->perfRateBrine_[this->first_perf_index_[well_index]];
}
/// One rate pr well
double brineWellRate(const int w) const;
@@ -238,28 +258,28 @@ public:
return well_potentials_;
}
std::vector<double>& perfThroughput() {
return perf_water_throughput_;
double * perfThroughput(std::size_t well_index) {
return &perf_water_throughput_[this->first_perf_index_[well_index]];
}
const std::vector<double>& perfThroughput() const {
return perf_water_throughput_;
const double * perfThroughput(std::size_t well_index) const {
return &perf_water_throughput_[this->first_perf_index_[well_index]];
}
std::vector<double>& perfSkinPressure() {
return perf_skin_pressure_;
double * perfSkinPressure(std::size_t well_index) {
return &perf_skin_pressure_[this->first_perf_index_[well_index]];
}
const std::vector<double>& perfSkinPressure() const {
return perf_skin_pressure_;
const double * perfSkinPressure(std::size_t well_index) const {
return &perf_skin_pressure_[this->first_perf_index_[well_index]];
}
std::vector<double>& perfWaterVelocity() {
return perf_water_velocity_;
double * perfWaterVelocity(std::size_t well_index) {
return &perf_water_velocity_[this->first_perf_index_[well_index]];
}
const std::vector<double>& perfWaterVelocity() const {
return perf_water_velocity_;
const double * perfWaterVelocity(std::size_t well_index) const {
return &perf_water_velocity_[this->first_perf_index_[well_index]];
}
template<class Comm>