OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-13 09:51:57 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3184 from totto82/fixZeroInit

Browse Source 
Improve initialization of the well rates for trivial rates
This commit is contained in:
Tor Harald Sandve 2021-05-10 08:53:16 +02:00 committed by GitHub
commit c87c2666d1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 169 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -422,6 +422,7 @@ namespace Opm {
const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
if (event) {
try {
well->updateWellStateWithTarget(ebosSimulator_, this->wellState(), local_deferredLogger);
well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), local_deferredLogger);
well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
} catch (const std::exception& e) {

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

@ -521,6 +521,8 @@ namespace Opm
double scalingFactor(const int comp_idx) const;
std::vector<double> initialWellRateFrations(const Simulator& ebosSimulator, const std::vector<double>& potentials) const;
// whether a well is specified with a non-zero and valid VFP table number
bool isVFPActive(Opm::DeferredLogger& deferred_logger) const;

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

@ -1695,11 +1695,33 @@ namespace Opm
}
double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
well_state.bhp()[well_index] = bhp;
// if the total rates are negative or zero
// 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*np + p] = well_state.wellPotentials()[well_index*np + p];
}
}
break;
}
case Well::InjectorCMode::BHP:
{
well_state.bhp()[well_index] = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate += well_state.wellRates()[well_index*np + 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*np + p] = well_state.wellPotentials()[well_index*np + p];
}
}
break;
}
case Well::InjectorCMode::GRUP:
@ -1719,42 +1741,66 @@ namespace Opm
{
const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
const auto& controls = well.productionControls(summaryState);
switch (current) {
case Well::ProducerCMode::ORAT:
{
double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Oil] ];
if (current_rate == 0.0)
break;
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= controls.oil_rate/current_rate;
// 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*np + p] *= controls.oil_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
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*np + p] = - fractions[p] * controls.oil_rate/control_fraction;
}
}
}
break;
}
case Well::ProducerCMode::WRAT:
{
double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Water] ];
if (current_rate == 0.0)
break;
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= controls.water_rate/current_rate;
// 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*np + p] *= controls.water_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
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*np + p] = - fractions[p] * controls.water_rate/control_fraction;
}
}
}
break;
}
case Well::ProducerCMode::GRAT:
{
double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Gas] ];
if (current_rate == 0.0)
break;
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= controls.gas_rate/current_rate;
// 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*np + p] *= controls.gas_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
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*np + p] = - fractions[p] * controls.gas_rate/control_fraction;
}
}
}
break;
}
@ -1762,12 +1808,20 @@ namespace Opm
{
double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Water] ]
- well_state.wellRates()[ well_index*np + pu.phase_pos[Oil] ];
if (current_rate == 0.0)
break;
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= controls.liquid_rate/current_rate;
// 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*np + p] *= controls.liquid_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
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*np + p] = - fractions[p] * controls.liquid_rate / control_fraction;
}
}
}
break;
}
@ -1783,12 +1837,18 @@ namespace Opm
for (int p = 0; p<np; ++p) {
total_res_rate -= well_state.wellRates()[well_index*np + p] * convert_coeff[p];
}
if (total_res_rate == 0.0)
break;
if (controls.prediction_mode) {
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= controls.resv_rate/total_res_rate;
// 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*np + p] *= controls.resv_rate/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
}
}
} else {
std::vector<double> hrates(number_of_phases_,0.);
@ -1804,8 +1864,17 @@ namespace Opm
std::vector<double> hrates_resv(number_of_phases_,0.);
rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, pvtRegionIdx_, hrates, hrates_resv);
double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] *= target/total_res_rate;
// 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*np + p] *= target/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFrations(ebos_simulator, well_state.wellPotentials());
for (int p = 0; p<np; ++p) {
well_state.wellRates()[well_index*np + p] = - fractions[p] * target / convert_coeff[p];
}
}
}
@ -1814,6 +1883,18 @@ namespace Opm
case Well::ProducerCMode::BHP:
{
well_state.bhp()[well_index] = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate -= well_state.wellRates()[well_index*np + 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*np + p] = -well_state.wellPotentials()[well_index*np + p];
}
}
break;
}
case Well::ProducerCMode::THP:
@ -1824,6 +1905,17 @@ namespace Opm
}
double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
well_state.bhp()[well_index] = bhp;
// if the total rates are negative or zero
// 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*np + p] = -well_state.wellPotentials()[well_index*np + p];
}
}
break;
}
case Well::ProducerCMode::GRUP:
{
@ -1841,6 +1933,48 @@ namespace Opm
}
}
template<typename TypeTag>
std::vector<double>
WellInterface<TypeTag>::
initialWellRateFrations(const Simulator& ebosSimulator, const std::vector<double>& potentials) const
{
const int np = number_of_phases_;
std::vector<double> scaling_factor(np);
double total_potentials = 0.0;
for (int p = 0; p<np; ++p) {
total_potentials += potentials[index_of_well_*np + p];
}
if (total_potentials > 0) {
for (int p = 0; p<np; ++p) {
scaling_factor[p] = potentials[index_of_well_*np + p] / total_potentials;
}
return scaling_factor;
}
// if we don't have any potentials we weight it using the mobilites
// We only need approximation so we don't bother with the vapporized oil and dissolved gas
double total_tw = 0;
const int nperf = number_of_perforations_;
for (int perf = 0; perf < nperf; ++perf) {
total_tw += well_index_[perf];
}
for (int perf = 0; perf < nperf; ++perf) {
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
const double well_tw_fraction = well_index_[perf] / total_tw;
double total_mobility = 0.0;
for (int p = 0; p < np; ++p) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(p);
total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
}
for (int p = 0; p < np; ++p) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(p);
scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
}
}
return scaling_factor;
}
template<typename TypeTag>
bool