OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3336 from joakim-hove/wellpotential_wellcontainer

Browse Source 
Manage well potentials with WellContainer<>
This commit is contained in:
Joakim Hove 2021-06-06 08:33:26 +02:00 committed by GitHub
commit e83d4d0dbd
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 64 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -1619,7 +1619,7 @@ namespace Opm {
// potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
// and updated only if sucessfull. i.e. the potentials are zero for exceptions
for (int p = 0; p < np; ++p) {
this->wellState().wellPotentials()[well->indexOfWell() * np + p] = std::abs(potentials[p]);
this->wellState().wellPotentials(well->indexOfWell())[p] = std::abs(potentials[p]);
}
}
}
@ -2787,8 +2787,7 @@ namespace Opm {
wellPI(const int well_index) const
{
const auto& pu = this->phase_usage_;
const auto np = this->numPhases();
const auto* pi = &this->wellState().productivityIndex()[np*well_index + 0];
const auto& pi = this->wellState().productivityIndex(well_index);
const auto preferred = this->wells_ecl_[well_index].getPreferredPhase();
switch (preferred) { // Should really have LIQUID = OIL + WATER here too...

6
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -632,7 +632,7 @@ namespace Opm
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase){
well_state_copy.wellRates(well_copy.index_of_well_)[phase]
= sign * well_state_copy.wellPotentials()[well_copy.index_of_well_*np + phase];
= sign * well_state_copy.wellPotentials(well_copy.index_of_well_)[phase];
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
@ -992,8 +992,8 @@ namespace Opm
std::transform(src, src + np, dest, dest, std::plus<>{});
};
auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
auto* wellPI = well_state.productivityIndex(this->index_of_well_).data();
auto* connPI = well_state.connectionProductivityIndex(this->index_of_well_).data();
setToZero(wellPI);

4
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -2042,8 +2042,8 @@ namespace Opm
std::transform(src, src + np, dest, dest, std::plus<>{});
};
auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
auto* wellPI = well_state.productivityIndex(this->index_of_well_).data();
auto* connPI = well_state.connectionProductivityIndex(this->index_of_well_).data();
setToZero(wellPI);

5
opm/simulators/wells/WellGroupHelpers.cpp
View File

@ -1279,10 +1279,9 @@ namespace WellGroupHelpers
continue;
}
const auto wellrate_index = well_index * wellState.numPhases();
// add contribution from wells unconditionally
for (int phase = 0; phase < np; phase++) {
pot[phase] += wefac * wellState.wellPotentials()[wellrate_index + phase];
pot[phase] += wefac * wellState.wellPotentials(well_index)[phase];
}
}
@ -1332,7 +1331,7 @@ namespace WellGroupHelpers
// the well is found and owned
int well_index = it->second[0];
const auto wpot = wellState.wellPotentials().data() + well_index * wellState.numPhases();
const auto wpot = wellState.wellPotentials(well_index);
if (pu.phase_used[BlackoilPhases::Liquid] > 0)
oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]];

2
opm/simulators/wells/WellInterface.hpp
View File

@ -314,7 +314,7 @@ protected:
double scalingFactor(const int comp_idx) const;
std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const std::vector<double>& potentials) const;
std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const;
// check whether the well is operable under BHP limit with current reservoir condition
virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) =0;

3
opm/simulators/wells/WellInterfaceFluidSystem.cpp
View File

@ -685,10 +685,9 @@ updateWellTestStateEconomic(const WellState& well_state,
bool rate_limit_violated = false;
const auto& quantity_limit = econ_production_limits.quantityLimit();
const int np = number_of_phases_;
if (econ_production_limits.onAnyRateLimit()) {
if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN)
rate_limit_violated = checkRateEconLimits(econ_production_limits, &well_state.wellPotentials()[index_of_well_ * np], deferred_logger);
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellPotentials(index_of_well_).data(), deferred_logger);
else {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger);
}

26
opm/simulators/wells/WellInterface_impl.hpp
View File

@ -745,7 +745,7 @@ namespace Opm
double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = well_state.wellPotentials()[well_index*np + p];
well_state.wellRates(well_index)[p] = well_state.wellPotentials(well_index)[p];
}
}
break;
@ -762,7 +762,7 @@ namespace Opm
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = well_state.wellPotentials()[well_index*np + p];
well_state.wellRates(well_index)[p] = well_state.wellPotentials(well_index)[p];
}
}
break;
@ -795,7 +795,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= controls.oil_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
@ -815,7 +815,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= controls.water_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
@ -835,7 +835,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= controls.gas_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Gas]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
@ -858,7 +858,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= controls.liquid_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
@ -888,7 +888,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= controls.resv_rate/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
}
@ -914,7 +914,7 @@ namespace Opm
well_state.wellRates(well_index)[p] *= target/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state.wellPotentials());
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * target / convert_coeff[p];
}
@ -935,7 +935,7 @@ namespace Opm
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -well_state.wellPotentials()[well_index*np + p];
well_state.wellRates(well_index)[p] = -well_state.wellPotentials(well_index)[p];
}
}
break;
@ -955,7 +955,7 @@ namespace Opm
double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -well_state.wellPotentials()[well_index*np + p];
well_state.wellRates(well_index)[p] = -well_state.wellPotentials(well_index)[p];
}
}
break;
@ -979,18 +979,18 @@ namespace Opm
template<typename TypeTag>
std::vector<double>
WellInterface<TypeTag>::
initialWellRateFractions(const Simulator& ebosSimulator, const std::vector<double>& potentials) const
initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const
{
const int np = this->number_of_phases_;
std::vector<double> scaling_factor(np);
double total_potentials = 0.0;
for (int p = 0; p<np; ++p) {
total_potentials += potentials[this->index_of_well_*np + p];
total_potentials += well_state.wellPotentials(this->index_of_well_)[p];
}
if (total_potentials > 0) {
for (int p = 0; p<np; ++p) {
scaling_factor[p] = potentials[this->index_of_well_*np + p] / total_potentials;
scaling_factor[p] = well_state.wellPotentials(this->index_of_well_)[p] / total_potentials;
}
return scaling_factor;
}

62
opm/simulators/wells/WellState.cpp
View File

@ -57,6 +57,9 @@ void WellState::base_init(const std::vector<double>& cellPressures,
this->thp_.clear();
this->temperature_.clear();
this->segment_state.clear();
this->well_potentials_.clear();
this->productivity_index_.clear();
this->conn_productivity_index_.clear();
{
// const int nw = wells->number_of_wells;
const int nw = wells_ecl.size();
@ -103,7 +106,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
this->well_perf_data_.add(well.name(), well_perf_data);
this->parallel_well_info_.add(well.name(), well_info);
this->wellrates_.add(well.name(), std::vector<double>(np, 0));
this->well_potentials_.add(well.name(), std::vector<double>(np, 0));
const int num_perf_this_well = well_info->communication().sum(well_perf_data_[w].size());
this->segment_state.add(well.name(), SegmentState{});
this->perfpress_.add(well.name(), std::vector<double>(num_perf_this_well, -1e100));
@ -117,6 +120,8 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
this->perfRateBrine_.add(well.name(), std::vector<double>(num_perf_this_well, 0));
this->bhp_.add(well.name(), 0.0);
this->thp_.add(well.name(), 0.0);
this->productivity_index_.add(well.name(), std::vector<double>(np, 0));
this->conn_productivity_index_.add(well.name(), std::vector<double>(num_perf_this_well * np, 0));
if ( well.isInjector() )
this->temperature_.add(well.name(), well.injectionControls(summary_state).temperature);
else
@ -341,9 +346,6 @@ void WellState::init(const std::vector<double>& cellPressures,
}
}
productivity_index_.assign(nw * this->numPhases(), 0.0);
conn_productivity_index_.assign(nperf * this->numPhases(), 0.0);
well_potentials_.assign(nw * this->numPhases(), 0.0);
for (int w = 0; w < nw; ++w) {
switch (wells_ecl[w].getStatus()) {
@ -371,7 +373,6 @@ void WellState::init(const std::vector<double>& cellPressures,
continue;
}
const auto& wname = well.name();
auto it = prevState->wellMap().find(well.name());
if (it != end)
{
@ -403,9 +404,8 @@ void WellState::init(const std::vector<double>& cellPressures,
wellReservoirRates(w) = prevState->wellReservoirRates(oldIndex);
// Well potentials
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
wellPotentials()[ idx ] = prevState->wellPotentials()[ oldidx ];
for (int p=0; p < np; p++) {
this->wellPotentials(newIndex)[p] = prevState->wellPotentials(oldIndex)[p];
}
// perfPhaseRates
@ -481,15 +481,8 @@ void WellState::init(const std::vector<double>& cellPressures,
}
}
// Productivity index.
{
auto* thisWellPI = &this ->productivityIndex()[newIndex*np + 0];
const auto* thatWellPI = &prevState->productivityIndex()[oldIndex*np + 0];
for (int p = 0; p < np; ++p) {
thisWellPI[p] = thatWellPI[p];
}
}
// Productivity index.
this->productivity_index_.copy_welldata( prevState->productivity_index_, wname );
}
// If in the new step, there is no THP related
@ -644,28 +637,26 @@ WellState::report(const int* globalCellIdxMap,
}
auto& well = res.at(wt.first);
const int well_rate_index = w * pu.num_phases;
const auto& reservoir_rates = this->well_reservoir_rates_[w];
const auto& well_potentials = this->well_potentials_[w];
const auto& wpi = this->productivity_index_[w];
if (pu.phase_used[Water]) {
const auto i = well_rate_index + pu.phase_pos[Water];
well.rates.set(rt::reservoir_water, reservoir_rates[pu.phase_pos[Water]]);
well.rates.set(rt::productivity_index_water, this->productivity_index_[i]);
well.rates.set(rt::well_potential_water, this->well_potentials_[i]);
well.rates.set(rt::productivity_index_water, wpi[pu.phase_pos[Water]]);
well.rates.set(rt::well_potential_water, well_potentials[pu.phase_pos[Water]]);
}
if (pu.phase_used[Oil]) {
const auto i = well_rate_index + pu.phase_pos[Oil];
well.rates.set(rt::reservoir_oil, reservoir_rates[pu.phase_pos[Oil]]);
well.rates.set(rt::productivity_index_oil, this->productivity_index_[i]);
well.rates.set(rt::well_potential_oil, this->well_potentials_[i]);
well.rates.set(rt::productivity_index_oil, wpi[pu.phase_pos[Oil]]);
well.rates.set(rt::well_potential_oil, well_potentials[pu.phase_pos[Oil]]);
}
if (pu.phase_used[Gas]) {
const auto i = well_rate_index + pu.phase_pos[Gas];
well.rates.set(rt::reservoir_gas, reservoir_rates[pu.phase_pos[Gas]]);
well.rates.set(rt::productivity_index_gas, this->productivity_index_[i]);
well.rates.set(rt::well_potential_gas, this->well_potentials_[i]);
well.rates.set(rt::productivity_index_gas, wpi[pu.phase_pos[Gas]]);
well.rates.set(rt::well_potential_gas, well_potentials[pu.phase_pos[Gas]]);
}
if (pu.has_solvent || pu.has_zFraction) {
@ -752,14 +743,12 @@ void WellState::reportConnections(data::Well& well,
pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas;
}
for( auto& comp : well.connections) {
const auto connPhaseOffset = np * (wt.second[1] + local_comp_index);
const auto& rates = &this->perfPhaseRates(well_index)[np*local_comp_index];
const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset;
const auto * rates = &this->perfPhaseRates(well_index)[np*local_comp_index];
const auto& connPI = this->connectionProductivityIndex(well_index);
for( int i = 0; i < np; ++i ) {
comp.rates.set( phs[ i ], rates[i] );
comp.rates.set( pi [ i ], *(connPI + i) );
comp.rates.set( pi [ i ], connPI[i] );
}
if ( pu.has_polymer ) {
const auto& perf_polymer_rate = this->perfRatePolymer(well_index);
@ -794,7 +783,6 @@ void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
// what we do here, is to set the segment rates and perforation rates
for (int w = 0; w < nw; ++w) {
const auto& well_ecl = wells_ecl[w];
const auto& wname = wells_ecl[w].name();
if ( well_ecl.isMultiSegment() ) {
const WellSegments& segment_set = well_ecl.getSegments();
@ -966,17 +954,15 @@ void WellState::shutWell(int well_index)
this->wellrates_[well_index].assign(np, 0);
auto& resv = this->well_reservoir_rates_[well_index];
auto* wpi = &this->productivity_index_[np*well_index + 0];
auto& wpi = this->productivity_index_[well_index];
for (int p = 0; p < np; ++p) {
resv[p] = 0.0;
wpi[p] = 0.0;
}
const auto first = this->first_perf_index_[well_index]*np;
const auto last = first + this->numPerf(well_index)*np;
std::fill(this->conn_productivity_index_.begin() + first,
this->conn_productivity_index_.begin() + last, 0.0);
auto& connpi = this->conn_productivity_index_[well_index];
connpi.assign(connpi.size(), 0);
}
void WellState::updateStatus(int well_index, Well::Status status)

31
opm/simulators/wells/WellState.hpp
View File

@ -235,28 +235,29 @@ public:
return this->segment_state[wname];
}
std::vector<double>& productivityIndex() {
return productivity_index_;
std::vector<double>& productivityIndex(std::size_t well_index) {
return this->productivity_index_[well_index];
}
const std::vector<double>& productivityIndex() const {
return productivity_index_;
const std::vector<double>& productivityIndex(std::size_t well_index) const {
return this->productivity_index_[well_index];
}
std::vector<double>& connectionProductivityIndex() {
return this->conn_productivity_index_;
std::vector<double>& connectionProductivityIndex(std::size_t well_index) {
return this->conn_productivity_index_[well_index];
}
const std::vector<double>& connectionProductivityIndex() const {
return this->conn_productivity_index_;
const std::vector<double>& connectionProductivityIndex(std::size_t well_index) const {
return this->conn_productivity_index_[well_index];
}
std::vector<double>& wellPotentials() {
return well_potentials_;
std::vector<double>& wellPotentials(std::size_t well_index) {
return this->well_potentials_[well_index];
}
const std::vector<double>& wellPotentials() const {
return well_potentials_;
const std::vector<double>& wellPotentials(std::size_t well_index) const {
return this->well_potentials_[well_index];
}
std::vector<double>& perfThroughput(std::size_t well_index) {
@ -492,13 +493,13 @@ private:
WellContainer<SegmentState> segment_state;
// Productivity Index
std::vector<double> productivity_index_;
WellContainer<std::vector<double>> productivity_index_;
// Connection-level Productivity Index
std::vector<double> conn_productivity_index_;
WellContainer<std::vector<double>> conn_productivity_index_;
// Well potentials
std::vector<double> well_potentials_;
WellContainer<std::vector<double>> well_potentials_;
data::Segment