Merge pull request #3277 from joakim-hove/segment-index

Segment index
This commit is contained in:
Joakim Hove 2021-05-20 10:01:07 +02:00 committed by GitHub
commit e5ab36e0f1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 42 additions and 35 deletions

View File

@ -2013,21 +2013,23 @@ namespace Opm {
const WellSegments& segment_set = well_ecl.getSegments(); const WellSegments& segment_set = well_ecl.getSegments();
const int top_segment_index = well_state.topSegmentIndex(well_index); const int top_segment_index = well_state.topSegmentIndex(well_index);
const auto& segments = rst_well.segments; const auto& rst_segments = rst_well.segments;
// \Note: eventually we need to hanlde the situations that some segments are shut // \Note: eventually we need to hanlde the situations that some segments are shut
assert(0u + segment_set.size() == segments.size()); assert(0u + segment_set.size() == rst_segments.size());
for (const auto& segment : segments) { auto * segment_pressure = &well_state.segPress()[top_segment_index];
const int segment_index = segment_set.segmentNumberToIndex(segment.first); auto * segment_rates = &well_state.segRates()[top_segment_index*np];
for (const auto& rst_segment : rst_segments) {
const int segment_index = segment_set.segmentNumberToIndex(rst_segment.first);
// recovering segment rates and pressure from the restart values // recovering segment rates and pressure from the restart values
const auto pres_idx = data::SegmentPressures::Value::Pressure; const auto pres_idx = data::SegmentPressures::Value::Pressure;
well_state.segPress()[top_segment_index + segment_index] = segment.second.pressures[pres_idx]; segment_pressure[segment_index] = rst_segment.second.pressures[pres_idx];
const auto& segment_rates = segment.second.rates; const auto& rst_segment_rates = rst_segment.second.rates;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
well_state.segRates()[(top_segment_index + segment_index) * np + p] = segment_rates.get(phs[p]); segment_rates[segment_index * np + p] = rst_segment_rates.get(phs[p]);
} }
} }
} }

View File

@ -279,14 +279,14 @@ namespace Opm
scaleSegmentRatesWithWellRates(WellState& well_state) const scaleSegmentRatesWithWellRates(WellState& well_state) const
{ {
const int top_segment_index = well_state.topSegmentIndex(index_of_well_); const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
auto * segment_rates = &well_state.segRates()[top_segment_index*this->number_of_phases_];
for (int phase = 0; phase < number_of_phases_; ++phase) { for (int phase = 0; phase < number_of_phases_; ++phase) {
const double well_phase_rate = well_state.wellRates()[number_of_phases_*index_of_well_ + phase]; const double well_phase_rate = well_state.wellRates()[number_of_phases_*index_of_well_ + phase];
const double unscaled_top_seg_rate = well_state.segRates()[number_of_phases_*top_segment_index + phase]; const double unscaled_top_seg_rate = segment_rates[phase];
if (std::abs(unscaled_top_seg_rate) > 1e-12) if (std::abs(unscaled_top_seg_rate) > 1e-12)
{ {
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
const int seg_index = top_segment_index + seg; segment_rates[this->number_of_phases_*seg + phase] *= well_phase_rate/unscaled_top_seg_rate;
well_state.segRates()[number_of_phases_*seg_index + phase] *= well_phase_rate/unscaled_top_seg_rate;
} }
} else { } else {
// for newly opened wells, the unscaled rate top segment rate is zero // for newly opened wells, the unscaled rate top segment rate is zero
@ -302,11 +302,10 @@ namespace Opm
perforation_rates[number_of_phases_ * perf + phase] = well_index_[perf] * perf_phaserate_scaled; perforation_rates[number_of_phases_ * perf + phase] = well_index_[perf] * perf_phaserate_scaled;
} }
std::vector<double> segment_rates; std::vector<double> rates;
WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_, WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_,
0, segment_rates); 0, rates);
std::copy(segment_rates.begin(), segment_rates.end(), std::copy(rates.begin(), rates.end(), segment_rates);
well_state.segRates().begin() + number_of_phases_ * top_segment_index );
} }
} }
} }
@ -319,10 +318,10 @@ namespace Opm
//scale segment pressures //scale segment pressures
const int top_segment_index = well_state.topSegmentIndex(index_of_well_); const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
const double bhp = well_state.bhp(index_of_well_); const double bhp = well_state.bhp(index_of_well_);
const double unscaled_top_seg_pressure = well_state.segPress()[top_segment_index]; auto * segment_pressure = &well_state.segPress()[top_segment_index];
const double unscaled_top_seg_pressure = segment_pressure[0];
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
const int seg_index = top_segment_index + seg; segment_pressure[seg] *= bhp/unscaled_top_seg_pressure;
well_state.segPress()[seg_index] *= bhp/unscaled_top_seg_pressure;
} }
} }
@ -724,30 +723,30 @@ namespace Opm
// the index of the top segment in the WellState // the index of the top segment in the WellState
const int top_segment_index = well_state.topSegmentIndex(index_of_well_); const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
const std::vector<double>& segment_rates = well_state.segRates(); const auto * segment_rates = &well_state.segRates()[top_segment_index * this->number_of_phases_];
const auto * segment_pressure = &well_state.segPress()[top_segment_index];
const PhaseUsage& pu = phaseUsage(); const PhaseUsage& pu = phaseUsage();
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
// calculate the total rate for each segment // calculate the total rate for each segment
double total_seg_rate = 0.0; double total_seg_rate = 0.0;
const int seg_index = top_segment_index + seg;
// the segment pressure // the segment pressure
primary_variables_[seg][SPres] = well_state.segPress()[seg_index]; primary_variables_[seg][SPres] = segment_pressure[seg];
// TODO: under what kind of circustances, the following will be wrong? // TODO: under what kind of circustances, the following will be wrong?
// the definition of g makes the gas phase is always the last phase // the definition of g makes the gas phase is always the last phase
for (int p = 0; p < number_of_phases_; p++) { for (int p = 0; p < number_of_phases_; p++) {
total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p]; total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg + p];
} }
primary_variables_[seg][GTotal] = total_seg_rate; primary_variables_[seg][GTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) { if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) { if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate; primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg + water_pos] / total_seg_rate;
} }
if (has_gfrac_variable) { if (has_gfrac_variable) {
const int gas_pos = pu.phase_pos[Gas]; const int gas_pos = pu.phase_pos[Gas];
primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate; primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg + gas_pos] / total_seg_rate;
} }
} else { // total_seg_rate == 0 } else { // total_seg_rate == 0
if (this->isInjector()) { if (this->isInjector()) {
@ -2341,6 +2340,8 @@ namespace Opm
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil]; const int oil_pos = pu.phase_pos[Oil];
auto * segment_rates = &well_state.segRates()[well_state.topSegmentIndex(this->index_of_well_) * this->number_of_phases_];
auto * segment_pressure = &well_state.segPress()[well_state.topSegmentIndex(this->index_of_well_)];
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
std::vector<double> fractions(number_of_phases_, 0.0); std::vector<double> fractions(number_of_phases_, 0.0);
fractions[oil_pos] = 1.0; fractions[oil_pos] = 1.0;
@ -2371,19 +2372,18 @@ namespace Opm
// calculate the phase rates based on the primary variables // calculate the phase rates based on the primary variables
const double g_total = primary_variables_[seg][GTotal]; const double g_total = primary_variables_[seg][GTotal];
const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
for (int p = 0; p < number_of_phases_; ++p) { for (int p = 0; p < number_of_phases_; ++p) {
const double phase_rate = g_total * fractions[p]; const double phase_rate = g_total * fractions[p];
well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate; segment_rates[seg*this->number_of_phases_ + p] = phase_rate;
if (seg == 0) { // top segment if (seg == 0) { // top segment
well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate; well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate;
} }
} }
// update the segment pressure // update the segment pressure
well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres]; segment_pressure[seg] = primary_variables_[seg][SPres];
if (seg == 0) { // top segment if (seg == 0) { // top segment
well_state.update_bhp(index_of_well_, well_state.segPress()[seg + top_segment_index]); well_state.update_bhp(index_of_well_, segment_pressure[seg]);
} }
} }
updateThp(well_state, deferred_logger); updateThp(well_state, deferred_logger);

View File

@ -685,21 +685,26 @@ void WellStateFullyImplicitBlackoil::initWellStateMSWell(const std::vector<Well>
} }
const int old_top_segment_index = prev_well_state->topSegmentIndex(old_index_well); const int old_top_segment_index = prev_well_state->topSegmentIndex(old_index_well);
const int new_top_segmnet_index = topSegmentIndex(new_index_well); const int new_top_segment_index = topSegmentIndex(new_index_well);
int number_of_segment = 0; int number_of_segment = 0;
// if it is the last well in list // if it is the last well in list
if (new_index_well == int(top_segment_index_.size()) - 1) { if (new_index_well == int(top_segment_index_.size()) - 1) {
number_of_segment = nseg_ - new_top_segmnet_index; number_of_segment = nseg_ - new_top_segment_index;
} else { } else {
number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index; number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segment_index;
} }
for (int i = 0; i < number_of_segment * np; ++i) { auto * segment_rates = &this->seg_rates_[new_top_segment_index*np];
seg_rates_[new_top_segmnet_index * np + i] = prev_well_state->segRates()[old_top_segment_index * np + i]; auto * segment_pressure = &this->seg_press_[new_top_segment_index];
}
for (int i = 0; i < number_of_segment; ++i) { const auto * prev_segment_rates = &prev_well_state->segRates()[old_top_segment_index*np];
seg_press_[new_top_segmnet_index + i] = prev_well_state->segPress()[old_top_segment_index + i]; const auto * prev_segment_pressure = &prev_well_state->segPress()[new_top_segment_index];
for (int seg=0; seg < number_of_segment; ++seg) {
for (int p = 0; p < np; ++p)
segment_rates[seg*np + p] = prev_segment_rates[seg*np + p];
segment_pressure[seg] = prev_segment_pressure[seg];
} }
} }
} }