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

Use SingleWellState to manage well surface rates

Browse Source
This commit is contained in:
Joakim Hove 2021-08-24 11:49:03 +02:00
parent 8937c9cba9
commit 755de65eb4
16 changed files with 136 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ebos/ecltracermodel.hh
View File

@ -434,7 +434,7 @@ protected:
// occasionally significant different from the sum over connections (as calculated above). Only observed
// for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations.
std::size_t well_index = simulator_.problem().wellModel().wellState().wellIndex(well.name());
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().wellRates(well_index)[tr.phaseIdx_];
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
rateWellTotal = official_well_rate_total;

6
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@ -186,7 +186,7 @@ loadRestartData(const data::Wells& rst_wells,
for( size_t i = 0; i < phs.size(); ++i ) {
assert( rst_well.rates.has( phs[ i ] ) );
well_state.wellRates(well_index)[i] = rst_well.rates.get(phs[i]);
ws.surface_rates[i] = rst_well.rates.get(phs[i]);
}
auto& perf_data = well_state.perfData(well_index);
@ -1173,7 +1173,7 @@ getGuideRateValues(const Well& well) const
auto grval = data::GuideRateValue{};
const auto& wname = well.name();
if (!this->wellState().hasWellRates(wname)) {
if (!this->wellState().has(wname)) {
// No flow rates for 'wname' -- might be before well comes
// online (e.g., for the initial condition before simulation
// starts).
@ -1333,7 +1333,7 @@ calculateAllGroupGuiderates(const int reportStepIdx) const
// group tree (FIELD group).
for (const auto& wname : schedule_.wellNames(reportStepIdx)) {
if (! (this->wellState().hasWellRates(wname) &&
if (! (this->wellState().has(wname) &&
this->guideRate_.has(wname)))
{
continue;

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

@ -300,7 +300,8 @@ GasLiftGroupInfo::
getProducerWellRates_(int well_index)
{
const auto& pu = this->phase_usage_;
const auto& wrate = this->well_state_.wellRates(well_index);
const auto& ws= this->well_state_.well(well_index);
const auto& wrate = ws.surface_rates;
const auto oil_rate = pu.phase_used[Oil]
? -wrate[pu.phase_pos[Oil]]

7
opm/simulators/wells/GasLiftStage2.cpp
View File

@ -404,11 +404,10 @@ GasLiftStage2::
getStdWellRates_(const WellInterfaceGeneric &well)
{
const int well_index = well.indexOfWell();
const auto& ws = this->well_state_.well(well_index);
const auto& pu = well.phaseUsage();
auto oil_rate =
-this->well_state_.wellRates(well_index)[pu.phase_pos[Oil]];
auto gas_rate =
-this->well_state_.wellRates(well_index)[pu.phase_pos[Gas]];
auto oil_rate = -ws.surface_rates[pu.phase_pos[Oil]];
auto gas_rate = -ws.surface_rates[pu.phase_pos[Gas]];
return {oil_rate, gas_rate};
}

16
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@ -1473,11 +1473,11 @@ updateThp(WellState& well_state,
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
auto& well = well_state.well(baseif_.indexOfWell());
well.thp = 0;
ws.thp = 0;
return;
}
@ -1486,18 +1486,16 @@ updateThp(WellState& well_state,
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Water ] ];
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Oil ] ];
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
auto& well = well_state.well(baseif_.indexOfWell());
const double bhp = well.bhp;
well.thp = this->calculateThpFromBhp(rates, bhp, rho, deferred_logger);
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
@ -1657,7 +1655,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
const double phase_rate = g_total * fractions[p];
segment_rates[seg*baseif_.numPhases() + p] = phase_rate;
if (seg == 0) { // top segment
well_state.wellRates(baseif_.indexOfWell())[p] = phase_rate;
ws.surface_rates[p] = phase_rate;
}
}

4
opm/simulators/wells/MultisegmentWellGeneric.cpp
View File

@ -105,7 +105,7 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
auto& segment_rates = segments.rates;
for (int phase = 0; phase < baseif_.numPhases(); ++phase) {
const double unscaled_top_seg_rate = segment_rates[phase];
const double well_phase_rate = well_state.wellRates(baseif_.indexOfWell())[phase];
const double well_phase_rate = ws.surface_rates[phase];
if (std::abs(unscaled_top_seg_rate) > 1e-12)
{
for (int seg = 0; seg < numberOfSegments(); ++seg) {
@ -120,7 +120,7 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
}
std::vector<double> perforation_rates(baseif_.numPhases() * baseif_.numPerfs(),0.0);
const double perf_phaserate_scaled = well_state.wellRates(baseif_.indexOfWell())[phase] / sumTw;
const double perf_phaserate_scaled = ws.surface_rates[phase] / sumTw;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
perforation_rates[baseif_.numPhases()* perf + phase] = baseif_.wellIndex()[perf] * perf_phaserate_scaled;
}

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

@ -279,14 +279,14 @@ namespace Opm
double total_rate = 0.0;
for (int phase = 0; phase < np; ++phase){
total_rate += well_state.wellRates(index_of_well_)[phase];
total_rate += ws.surface_rates[phase];
}
// for pressure controlled wells the well rates are the potentials
// if the rates are trivial we are most probably looking at the newly
// opened well and we therefore make the affort of computing the potentials anyway.
if (std::abs(total_rate) > 0) {
for (int phase = 0; phase < np; ++phase){
well_potentials[phase] = well_state.wellRates(index_of_well_)[phase];
well_potentials[phase] = ws.surface_rates[phase];
}
return;
}
@ -367,7 +367,7 @@ namespace Opm
const int np = number_of_phases_;
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 * ws.well_potentials[phase];
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);

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

@ -26,6 +26,7 @@ SingleWellState::SingleWellState(bool is_producer, std::size_t num_phases, doubl
, temperature(temp)
, well_potentials(num_phases)
, productivity_index(num_phases)
, surface_rates(num_phases)
{}
@ -49,6 +50,8 @@ void SingleWellState::shut() {
this->bhp = 0;
this->thp = 0;
this->status = Well::Status::SHUT;
std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0);
std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
}
void SingleWellState::stop() {

1
opm/simulators/wells/SingleWellState.hpp
View File

@ -41,6 +41,7 @@ public:
double vaporized_oil_rate{0};
std::vector<double> well_potentials;
std::vector<double> productivity_index;
std::vector<double> surface_rates;
SegmentState segments;
Events events;
Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};

46
opm/simulators/wells/StandardWellEval.cpp
View File

@ -255,11 +255,11 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
const int well_index = baseif_.indexOfWell();
const int np = baseif_.numPhases();
const auto& pu = baseif_.phaseUsage();
const auto& well = well_state.well(well_index);
const auto& ws = well_state.well(well_index);
// the weighted total well rate
double total_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
total_well_rate += baseif_.scalingFactor(p) * well_state.wellRates(well_index)[p];
total_well_rate += baseif_.scalingFactor(p) * ws.surface_rates[p];
}
// Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
@ -267,13 +267,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
if (baseif_.isInjector()) {
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Water]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
break;
case InjectorType::GAS:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Gas]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
break;
case InjectorType::OIL:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Oil]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
break;
case InjectorType::MULTI:
// Not supported.
@ -287,10 +287,10 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
if (std::abs(total_well_rate) > 0.) {
if constexpr (has_wfrac_variable) {
primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * well_state.wellRates(well_index)[pu.phase_pos[Water]] / total_well_rate;
primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
}
if constexpr (has_gfrac_variable) {
primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates(well_index)[pu.phase_pos[Gas]]
primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]]
- (Indices::enableSolvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ;
}
if constexpr (Indices::enableSolvent) {
@ -342,7 +342,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
// BHP
primary_variables_[Bhp] = well.bhp;
primary_variables_[Bhp] = ws.bhp;
}
template<class FluidSystem, class Indices, class Scalar>
@ -556,11 +556,11 @@ updateThp(WellState& well_state,
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
auto& well = well_state.well(baseif_.indexOfWell());
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
well.thp = 0;
ws.thp = 0;
return;
}
@ -569,20 +569,18 @@ updateThp(WellState& well_state,
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Water ] ];
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Oil ] ];
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
const double bhp = well.bhp;
well.thp = this->calculateThpFromBhp(well_state,
ws.thp = this->calculateThpFromBhp(well_state,
rates,
bhp,
ws.bhp,
deferred_logger);
}
@ -655,29 +653,29 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
F[pu.phase_pos[Gas]] += F_solvent;
}
auto& well = well_state.well(baseif_.indexOfWell());
well.bhp = primary_variables_[Bhp];
auto& ws = well_state.well(baseif_.indexOfWell());
ws.bhp = primary_variables_[Bhp];
// calculate the phase rates based on the primary variables
// for producers, this is not a problem, while not sure for injectors here
if (baseif_.isProducer()) {
const double g_total = primary_variables_[WQTotal];
for (int p = 0; p < baseif_.numPhases(); ++p) {
well_state.wellRates(baseif_.indexOfWell())[p] = g_total * F[p];
ws.surface_rates[p] = g_total * F[p];
}
} else { // injectors
for (int p = 0; p < baseif_.numPhases(); ++p) {
well_state.wellRates(baseif_.indexOfWell())[p] = 0.0;
ws.surface_rates[p] = 0.0;
}
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Water]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal];
break;
case InjectorType::GAS:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
break;
case InjectorType::OIL:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
break;
case InjectorType::MULTI:
// Not supported.

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

@ -1252,10 +1252,10 @@ namespace Opm
const int gaspos = gasCompIdx + perf * num_components_;
if (oilPresent) {
const double oilrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
if (oilrate > 0) {
const double gasrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
double rv = 0.0;
if (gasrate > 0) {
rv = oilrate / gasrate;
@ -1278,9 +1278,9 @@ namespace Opm
const int oilpos = oilCompIdx + perf * num_components_;
if (gasPresent) {
rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
const double gasrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
if (gasrate > 0) {
const double oilrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Oil]]);
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
double rs = 0.0;
if (oilrate > 0) {
rs = gasrate / oilrate;
@ -1862,14 +1862,14 @@ namespace Opm
double total_rate = 0.0;
for (int phase = 0; phase < np; ++phase){
total_rate += well_state.wellRates(index_of_well_)[phase];
total_rate += ws.surface_rates[phase];
}
// for pressure controlled wells the well rates are the potentials
// if the rates are trivial we are most probably looking at the newly
// opened well and we therefore make the affort of computing the potentials anyway.
if (std::abs(total_rate) > 0) {
for (int phase = 0; phase < np; ++phase){
well_potentials[phase] = well_state.wellRates(index_of_well_)[phase];
well_potentials[phase] = ws.surface_rates[phase];
}
return;
}

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

@ -97,6 +97,7 @@ namespace {
continue;
double factor = wellEcl.getEfficiencyFactor();
const auto& ws = wellState.well(well_index);
if (res_rates) {
const auto& well_rates = wellState.wellReservoirRates(well_index);
if (injector)
@ -104,7 +105,7 @@ namespace {
else
rate -= factor * well_rates[phasePos];
} else {
const auto& well_rates = wellState.wellRates(well_index);
const auto& well_rates = ws.surface_rates;
if (injector)
rate += factor * well_rates[phasePos];
else
@ -414,16 +415,17 @@ namespace WellGroupHelpers
const double efficiency = wellTmp.getEfficiencyFactor();
// add contributino from wells not under group control
const auto& ws_nupcol = wellStateNupcol.well(well_index);
const auto& ws = wellState.well(well_index);
if (isInjector) {
if (ws.injection_cmode != Well::InjectorCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase] += wellStateNupcol.wellRates(well_index)[phase] * efficiency;
groupTargetReduction[phase] += ws_nupcol.surface_rates[phase] * efficiency;
}
} else {
if (ws.production_cmode != Well::ProducerCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase] -= wellStateNupcol.wellRates(well_index)[phase] * efficiency;
groupTargetReduction[phase] -= ws_nupcol.surface_rates[phase] * efficiency;
}
}
}
@ -487,16 +489,16 @@ namespace WellGroupHelpers
}
// scale rates
const auto& ws = wellState.well(well_index);
auto& ws = wellState.well(well_index);
if (isInjector) {
if (ws.injection_cmode == Well::InjectorCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
wellState.wellRates(well_index)[phase] *= scale;
ws.surface_rates[phase] *= scale;
}
} else {
if (ws.production_cmode == Well::ProducerCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
wellState.wellRates(well_index)[phase] *= scale;
ws.surface_rates[phase] *= scale;
}
}
}
@ -578,8 +580,9 @@ namespace WellGroupHelpers
// opm-common that production and injection rates are positive.
if (!wellTmp.isInjector())
sign = -1;
const auto& ws = wellStateNupcol.well(well_index);
for (int phase = 0; phase < np; ++phase) {
rates[phase] = sign * wellStateNupcol.wellRates(well_index)[phase];
rates[phase] = sign * ws.surface_rates[phase];
}
}
wellState.setCurrentWellRates(wellName, rates);

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

@ -67,9 +67,10 @@ calculateReservoirRates(WellState& well_state) const
const int fipreg = 0; // not considering the region for now
const int np = number_of_phases_;
auto& ws = well_state.well(this->index_of_well_);
std::vector<double> surface_rates(np, 0.0);
for (int p = 0; p < np; ++p) {
surface_rates[p] = well_state.wellRates(index_of_well_)[p];
surface_rates[p] = ws.surface_rates[p];
}
std::vector<double> voidage_rates(np, 0.0);
@ -101,26 +102,26 @@ activeProductionConstraint(const WellState& well_state,
}
if (controls.hasControl(Well::ProducerCMode::ORAT) && currentControl != Well::ProducerCMode::ORAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Liquid]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
if (controls.oil_rate < current_rate)
return Well::ProducerCMode::ORAT;
}
if (controls.hasControl(Well::ProducerCMode::WRAT) && currentControl != Well::ProducerCMode::WRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Aqua]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (controls.water_rate < current_rate)
return Well::ProducerCMode::WRAT;
}
if (controls.hasControl(Well::ProducerCMode::GRAT) && currentControl != Well::ProducerCMode::GRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Vapour]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]];
if (controls.gas_rate < current_rate)
return Well::ProducerCMode::GRAT;
}
if (controls.hasControl(Well::ProducerCMode::LRAT) && currentControl != Well::ProducerCMode::LRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Aqua]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (controls.liquid_rate < current_rate)
return Well::ProducerCMode::LRAT;
}
@ -203,17 +204,17 @@ activeInjectionConstraint(const WellState& well_state,
switch (injectorType) {
case InjectorType::WATER:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Aqua] ];
break;
}
case InjectorType::OIL:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Liquid] ];
break;
}
case InjectorType::GAS:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Vapour] ];
break;
}
default:
@ -317,6 +318,7 @@ checkGroupConstraintsInj(const Group& group,
std::vector<double> resv_coeff(phaseUsage().num_phases, 1.0);
rateConverter_.calcInjCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(this->index_of_well_);
// Call check for the well's injection phase.
return WellGroupHelpers::checkGroupConstraintsInj(name(),
well_ecl_.groupName(),
@ -325,7 +327,7 @@ checkGroupConstraintsInj(const Group& group,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
ws.surface_rates.data(),
injectionPhase,
phaseUsage(),
efficiencyFactor,
@ -350,6 +352,7 @@ checkGroupConstraintsProd(const Group& group,
std::vector<double> resv_coeff(this->phaseUsage().num_phases, 1.0);
rateConverter_.calcCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(this->index_of_well_);
return WellGroupHelpers::checkGroupConstraintsProd(name(),
well_ecl_.groupName(),
group,
@ -357,7 +360,7 @@ checkGroupConstraintsProd(const Group& group,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
ws.surface_rates.data(),
phaseUsage(),
efficiencyFactor,
schedule,
@ -400,7 +403,7 @@ checkGroupConstraints(WellState& well_state,
ws.injection_cmode = Well::InjectorCMode::GRUP;
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
ws.surface_rates[p] *= group_constraint.second;
}
}
return group_constraint.first;
@ -428,7 +431,7 @@ checkGroupConstraints(WellState& well_state,
ws.production_cmode = Well::ProducerCMode::GRUP;
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
ws.surface_rates[p] *= group_constraint.second;
}
}
return group_constraint.first;
@ -698,7 +701,7 @@ updateWellTestStateEconomic(const WellState& well_state,
if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN)
rate_limit_violated = checkRateEconLimits(econ_production_limits, ws.well_potentials.data(), deferred_logger);
else {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger);
rate_limit_violated = checkRateEconLimits(econ_production_limits, ws.surface_rates.data(), deferred_logger);
}
}
@ -899,9 +902,9 @@ checkMaxRatioLimitWell(const WellState& well_state,
const int np = number_of_phases_;
std::vector<double> well_rates(np, 0.0);
const auto& ws = well_state.well(this->index_of_well_);
for (int p = 0; p < np; ++p) {
well_rates[p] = well_state.wellRates(index_of_well_)[p];
well_rates[p] = ws.surface_rates[p];
}
const double well_ratio = ratioFunc(well_rates, phaseUsage());
@ -1099,7 +1102,8 @@ getGroupProductionTargetRate(const Group& group,
}
// Avoid negative target rates coming from too large local reductions.
const double target_rate = std::max(0.0, target / efficiencyFactor);
const auto& rates = well_state.wellRates(index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& rates = ws.surface_rates;
const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
double scale = 1.0;
if (current_rate > 1e-14)

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

@ -610,7 +610,7 @@ namespace Opm
if (this->wellIsStopped()) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = 0.0;
ws.surface_rates[p] = 0;
}
ws.thp = 0;
return;
@ -647,7 +647,7 @@ namespace Opm
switch(current) {
case Well::InjectorCMode::RATE:
{
well_state.wellRates(well_index)[phasePos] = controls.surface_rate;
ws.surface_rates[phasePos] = controls.surface_rate;
break;
}
@ -656,16 +656,13 @@ namespace Opm
std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
const double coeff = convert_coeff[phasePos];
well_state.wellRates(well_index)[phasePos] = controls.reservoir_rate/coeff;
ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
break;
}
case Well::InjectorCMode::THP:
{
std::vector<double> rates(3, 0.0);
for (int p = 0; p<np; ++p) {
rates[p] = well_state.wellRates(well_index)[p];
}
auto rates = ws.surface_rates;
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
ws.bhp = bhp;
@ -673,11 +670,9 @@ namespace Opm
// we try to provide a better intial well rate
// using the well potentials
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] = ws.well_potentials[p];
}
}
if (total_rate <= 0.0)
ws.surface_rates = ws.well_potentials;
break;
}
case Well::InjectorCMode::BHP:
@ -685,16 +680,14 @@ namespace Opm
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate += well_state.wellRates(well_index)[p];
total_rate += ws.surface_rates[p];
}
// if the total rates are negative or zero
// we try to provide a better intial well rate
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = ws.well_potentials[p];
}
}
if (total_rate <= 0.0)
ws.surface_rates = ws.well_potentials;
break;
}
case Well::InjectorCMode::GRUP:
@ -712,7 +705,7 @@ namespace Opm
efficiencyFactor,
deferred_logger);
if (target)
well_state.wellRates(well_index)[phasePos] = *target;
ws.surface_rates[phasePos] = *target;
break;
}
case Well::InjectorCMode::CMODE_UNDEFINED:
@ -730,19 +723,19 @@ namespace Opm
switch (current) {
case Well::ProducerCMode::ORAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Oil] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.oil_rate/current_rate;
ws.surface_rates[p] *= controls.oil_rate/current_rate;
}
} else {
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) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.oil_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
}
}
}
@ -750,19 +743,19 @@ namespace Opm
}
case Well::ProducerCMode::WRAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Water] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.water_rate/current_rate;
ws.surface_rates[p] *= controls.water_rate/current_rate;
}
} else {
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) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.water_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
}
}
}
@ -770,19 +763,19 @@ namespace Opm
}
case Well::ProducerCMode::GRAT:
{
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[Gas] ];
double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.gas_rate/current_rate;
ws.surface_rates[p] *= controls.gas_rate/current_rate;
}
} else {
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) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.gas_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
}
}
}
@ -792,20 +785,20 @@ namespace Opm
}
case Well::ProducerCMode::LRAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Water] ]
- well_state.wellRates(well_index)[ pu.phase_pos[Oil] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
- ws.surface_rates[ pu.phase_pos[Oil] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.liquid_rate/current_rate;
ws.surface_rates[p] *= controls.liquid_rate/current_rate;
}
} else {
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) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.liquid_rate / control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
}
}
}
@ -821,19 +814,19 @@ namespace Opm
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
double total_res_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_res_rate -= well_state.wellRates(well_index)[p] * convert_coeff[p];
total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
}
if (controls.prediction_mode) {
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (total_res_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.resv_rate/total_res_rate;
ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
}
} else {
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];
ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
}
}
} else {
@ -854,12 +847,12 @@ namespace Opm
// but use either the potentials or the mobility ratio to initial the well rates
if (total_res_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= target/total_res_rate;
ws.surface_rates[p] *= target/total_res_rate;
}
} else {
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];
ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
}
}
@ -871,24 +864,21 @@ namespace Opm
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate -= well_state.wellRates(well_index)[p];
total_rate -= ws.surface_rates[p];
}
// if the total rates are negative or zero
// we try to provide a better intial well rate
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -ws.well_potentials[p];
ws.surface_rates[p] = -ws.well_potentials[p];
}
}
break;
}
case Well::ProducerCMode::THP:
{
std::vector<double> rates(3, 0.0);
for (int p = 0; p<np; ++p) {
rates[p] = well_state.wellRates(well_index)[p];
}
auto rates = ws.surface_rates;
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
ws.bhp = bhp;
@ -898,7 +888,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] = -ws.well_potentials[p];
ws.surface_rates[p] = -ws.well_potentials[p];
}
}
break;
@ -918,7 +908,7 @@ namespace Opm
// we don't want to scale with zero and get zero rates.
if (scale > 0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= scale;
ws.surface_rates[p] *= scale;
}
}
break;
@ -989,9 +979,10 @@ namespace Opm
{
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
auto& ws = well_state.well(this->index_of_well_);
int nonzero_rate_index = -1;
for (int p = 0; p < this->number_of_phases_; ++p) {
if (well_state.wellRates(this->index_of_well_)[p] != 0.0) {
if (ws.surface_rates[p] != 0.0) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {
@ -1009,12 +1000,12 @@ namespace Opm
std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
const double initial_nonzero_rate = well_state.wellRates(this->index_of_well_)[nonzero_rate_index];
const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
const int comp_idx_nz = this->flowPhaseToEbosCompIdx(nonzero_rate_index);
for (int p = 0; p < this->number_of_phases_; ++p) {
if (p != nonzero_rate_index) {
const int comp_idx = this->flowPhaseToEbosCompIdx(p);
double& rate = well_state.wellRates(this->index_of_well_)[p];
double& rate = ws.surface_rates[p];
rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
}
}

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

@ -42,7 +42,6 @@ void WellState::base_init(const std::vector<double>& cellPressures,
this->wellMap_.clear();
this->perfdata.clear();
this->parallel_well_info_.clear();
this->wellrates_.clear();
this->wells_.clear();
{
// const int nw = wells->number_of_wells;
@ -53,7 +52,7 @@ void WellState::base_init(const std::vector<double>& cellPressures,
const Well& well = wells_ecl[w];
// Initialize bhp(), thp(), wellRates(), temperature().
initSingleWell(cellPressures, w, well, well_perf_data[w], parallel_well_info[w], summary_state);
initSingleWell(cellPressures, well, well_perf_data[w], parallel_well_info[w], summary_state);
// Setup wellname -> well index mapping.
const int num_perf_this_well = well_perf_data[w].size();
@ -73,11 +72,10 @@ void WellState::base_init(const std::vector<double>& cellPressures,
void WellState::initSingleWell(const std::vector<double>& cellPressures,
const int w,
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,
const SummaryState& summary_state)
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,
const SummaryState& summary_state)
{
assert(well.isInjector() || well.isProducer());
@ -89,7 +87,6 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), static_cast<std::size_t>(np), temp});
this->parallel_well_info_.add(well.name(), well_info);
this->wellrates_.add(well.name(), std::vector<double>(np, 0));
const int num_perf_this_well = well_info->communication().sum(well_perf_data.size());
this->perfdata.add(well.name(), PerfData{static_cast<std::size_t>(num_perf_this_well), well.isInjector(), this->phase_usage_});
@ -142,7 +139,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
// (producer) or RATE (injector).
// Otherwise, we cannot set the correct
// value here and initialize to zero rate.
auto& rates = this->wellrates_[w];
auto& rates = ws.surface_rates;
if (well.isInjector()) {
if (inj_controls.cmode == Well::InjectorCMode::RATE) {
switch (inj_controls.injector_type) {
@ -267,6 +264,7 @@ void WellState::init(const std::vector<double>& cellPressures,
const int global_num_perf_this_well = ecl_well.getConnections().num_open();
auto& perf_data = this->perfData(w);
const auto& perf_input = well_perf_data[w];
const auto& ws = this->well(w);
for (int perf = 0; perf < num_perf_this_well; ++perf) {
perf_data.cell_index[perf] = perf_input[perf].cell_index;
@ -276,7 +274,7 @@ void WellState::init(const std::vector<double>& cellPressures,
if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
for (int p = 0; p < this->numPhases(); ++p) {
perf_data.phase_rates[this->numPhases()*perf + p] = wellRates(w)[p] / double(global_num_perf_this_well);
perf_data.phase_rates[this->numPhases()*perf + p] = ws.surface_rates[p] / double(global_num_perf_this_well);
}
}
perf_data.pressure[perf] = cellPressures[well_perf_data[w][perf].cell_index];
@ -348,7 +346,7 @@ void WellState::init(const std::vector<double>& cellPressures,
new_well.production_cmode = prev_well.production_cmode;
}
wellRates(w) = prevState->wellRates(oldIndex);
new_well.surface_rates = prev_well.surface_rates;
wellReservoirRates(w) = prevState->wellReservoirRates(oldIndex);
new_well.well_potentials = prev_well.well_potentials;
@ -380,7 +378,7 @@ void WellState::init(const std::vector<double>& cellPressures,
auto& target_rates = perf_data.phase_rates;
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);
target_rates[perf_index*np + p] = new_well.surface_rates[p] / double(global_num_perf_this_well);
}
}
}
@ -477,7 +475,7 @@ WellState::report(const int* globalCellIdxMap,
const auto& reservoir_rates = this->well_reservoir_rates_[well_index];
const auto& well_potentials = ws.well_potentials;
const auto& wpi = ws.productivity_index;
const auto& wv = this->wellRates(well_index);
const auto& wv = ws.surface_rates;
data::Well well;
well.bhp = ws.bhp;
@ -814,14 +812,10 @@ void WellState::shutWell(int well_index)
auto& ws = this->well(well_index);
ws.shut();
const int np = numPhases();
this->wellrates_[well_index].assign(np, 0);
auto& resv = this->well_reservoir_rates_[well_index];
auto& wpi = ws.productivity_index;
for (int p = 0; p < np; ++p) {
resv[p] = 0.0;
wpi[p] = 0.0;
}
auto& perf_data = this->perfData(well_index);

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

@ -255,8 +255,8 @@ public:
}
/// One rate per well and phase.
std::vector<double>& wellRates(std::size_t well_index) { return wellrates_[well_index]; }
const std::vector<double>& wellRates(std::size_t well_index) const { return wellrates_[well_index]; }
std::vector<double>& wellRates(std::size_t well_index) { return this->wells_[well_index].surface_rates; }
const std::vector<double>& wellRates(std::size_t well_index) const { return this->wells_[well_index].surface_rates; }
std::size_t numPerf(std::size_t well_index) const { return this->perfdata[well_index].size(); }
@ -312,7 +312,6 @@ private:
WellContainer<SingleWellState> wells_;
WellContainer<const ParallelWellInfo*> parallel_well_info_;
WellContainer<std::vector<double>> wellrates_;
WellContainer<PerfData> perfdata;
// The well_rates variable is defined for all wells on all processors. The
@ -358,7 +357,6 @@ private:
const SummaryState& summary_state);
void initSingleWell(const std::vector<double>& cellPressures,
const int w,
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,