Let the perforation index run over the indices for one well

This commit is contained in:
Joakim Hove 2021-05-11 19:35:26 +02:00
parent 846809ec29
commit b5580f39fa
4 changed files with 46 additions and 32 deletions

View File

@ -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.

View File

@ -584,6 +584,8 @@ 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];
for (int perf = 0; perf < number_of_perforations_; ++perf) {
// Calculate perforation quantities.
std::vector<EvalWell> 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<double> 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];

View File

@ -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()));

View File

@ -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<double> 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<double> perforation_rates(perf_rates, perf_rates + num_perf_this_well*np);
std::vector<double> 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