00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <opm/autodiff/StandardWells.hpp>
00023 #include <opm/autodiff/WellDensitySegmented.hpp>
00024
00025 #include <opm/autodiff/VFPInjProperties.hpp>
00026 #include <opm/autodiff/VFPProdProperties.hpp>
00027 #include <opm/autodiff/WellHelpers.hpp>
00028
00029
00030
00031
00032 namespace Opm
00033 {
00034
00035
00036 StandardWells::
00037 WellOps::WellOps(const Wells* wells)
00038 : w2p(),
00039 p2w(),
00040 well_cells()
00041 {
00042 if( wells )
00043 {
00044 w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
00045 p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
00046
00047 const int nw = wells->number_of_wells;
00048 const int* const wpos = wells->well_connpos;
00049
00050 typedef Eigen::Triplet<double> Tri;
00051
00052 std::vector<Tri> scatter, gather;
00053 scatter.reserve(wpos[nw]);
00054 gather .reserve(wpos[nw]);
00055
00056 for (int w = 0, i = 0; w < nw; ++w) {
00057 for (; i < wpos[ w + 1 ]; ++i) {
00058 scatter.push_back(Tri(i, w, 1.0));
00059 gather .push_back(Tri(w, i, 1.0));
00060 }
00061 }
00062
00063 w2p.setFromTriplets(scatter.begin(), scatter.end());
00064 p2w.setFromTriplets(gather .begin(), gather .end());
00065
00066 well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]);
00067 }
00068 }
00069
00070
00071
00072
00073
00074 StandardWells::StandardWells(const Wells* wells_arg, WellCollection* well_collection)
00075 : wells_active_(wells_arg!=nullptr)
00076 , wells_(wells_arg)
00077 , wops_(wells_arg)
00078 , well_collection_(well_collection)
00079 , well_perforation_efficiency_factors_(Vector::Ones(wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0))
00080 , fluid_(nullptr)
00081 , active_(nullptr)
00082 , phase_condition_(nullptr)
00083 , vfp_properties_(nullptr)
00084 , well_perforation_densities_(Vector())
00085 , well_perforation_pressure_diffs_(Vector())
00086 , store_well_perforation_fluxes_(false)
00087 {
00088 }
00089
00090
00091
00092
00093
00094 void
00095 StandardWells::init(const BlackoilPropsAdFromDeck* fluid_arg,
00096 const std::vector<bool>* active_arg,
00097 const std::vector<PhasePresence>* pc_arg,
00098 const VFPProperties* vfp_properties_arg,
00099 const double gravity_arg,
00100 const Vector& depth_arg)
00101 {
00102 fluid_ = fluid_arg;
00103 active_ = active_arg;
00104 phase_condition_ = pc_arg;
00105 vfp_properties_ = vfp_properties_arg;
00106 gravity_ = gravity_arg;
00107 perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
00108
00109 calculateEfficiencyFactors();
00110 }
00111
00112
00113
00114
00115
00116 const Wells& StandardWells::wells() const
00117 {
00118 assert(wells_ != 0);
00119 return *(wells_);
00120 }
00121
00122
00123 const Wells* StandardWells::wellsPointer() const
00124 {
00125 return wells_;
00126 }
00127
00128
00129
00130 bool StandardWells::wellsActive() const
00131 {
00132 return wells_active_;
00133 }
00134
00135
00136
00137
00138
00139 void StandardWells::setWellsActive(const bool wells_active)
00140 {
00141 wells_active_ = wells_active;
00142 }
00143
00144
00145
00146
00147
00148 bool StandardWells::localWellsActive() const
00149 {
00150 return wells_ ? (wells_->number_of_wells > 0 ) : false;
00151 }
00152
00153
00154
00155
00156
00157 int
00158 StandardWells::numWellVars() const
00159 {
00160 if ( !localWellsActive() )
00161 {
00162 return 0;
00163 }
00164
00165
00166 const int nw = wells().number_of_wells;
00167 return (numPhases() + 1) * nw;
00168 }
00169
00170
00171
00172
00173
00174 const StandardWells::WellOps&
00175 StandardWells::wellOps() const
00176 {
00177 return wops_;
00178 }
00179
00180
00181
00182
00183
00184 StandardWells::Vector& StandardWells::wellPerforationDensities()
00185 {
00186 return well_perforation_densities_;
00187 }
00188
00189
00190
00191
00192
00193 const StandardWells::Vector&
00194 StandardWells::wellPerforationDensities() const
00195 {
00196 return well_perforation_densities_;
00197 }
00198
00199
00200
00201
00202
00203 StandardWells::Vector&
00204 StandardWells::wellPerforationPressureDiffs()
00205 {
00206 return well_perforation_pressure_diffs_;
00207 }
00208
00209
00210
00211
00212
00213 const StandardWells::Vector&
00214 StandardWells::wellPerforationPressureDiffs() const
00215 {
00216 return well_perforation_pressure_diffs_;
00217 }
00218
00219
00220
00221
00222 template<class SolutionState, class WellState>
00223 void
00224 StandardWells::
00225 computePropertiesForWellConnectionPressures(const SolutionState& state,
00226 const WellState& xw,
00227 std::vector<double>& b_perf,
00228 std::vector<double>& rsmax_perf,
00229 std::vector<double>& rvmax_perf,
00230 std::vector<double>& surf_dens_perf)
00231 {
00232 const int nperf = wells().well_connpos[wells().number_of_wells];
00233 const int nw = wells().number_of_wells;
00234
00235
00236 const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf);
00237 Vector avg_press = perf_press*0;
00238 for (int w = 0; w < nw; ++w) {
00239 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
00240 const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
00241 const double p_avg = (perf_press[perf] + p_above)/2;
00242 avg_press[perf] = p_avg;
00243 }
00244 }
00245
00246 const std::vector<int>& well_cells = wellOps().well_cells;
00247
00248
00249 const ADB perf_temp = subset(state.temperature, well_cells);
00250
00251
00252
00253
00254 const ADB avg_press_ad = ADB::constant(avg_press);
00255 std::vector<PhasePresence> perf_cond(nperf);
00256
00257 for (int perf = 0; perf < nperf; ++perf) {
00258 perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
00259 }
00260 const PhaseUsage& pu = fluid_->phaseUsage();
00261 DataBlock b(nperf, pu.num_phases);
00262 if (pu.phase_used[BlackoilPhases::Aqua]) {
00263 const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
00264 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
00265 }
00266 assert((*active_)[Oil]);
00267 const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
00268 if (pu.phase_used[BlackoilPhases::Liquid]) {
00269 const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null();
00270 const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
00271 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
00272 }
00273 if (pu.phase_used[BlackoilPhases::Vapour]) {
00274 const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null();
00275 const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
00276 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
00277 }
00278 if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
00279 const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
00280 rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
00281
00282 const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
00283 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
00284 }
00285
00286
00287 b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
00288
00289
00290
00291
00292 Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
00293 for (int phase = 1; phase < pu.num_phases; ++phase) {
00294 rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
00295 }
00296 surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
00297
00298 }
00299
00300
00301
00302
00303
00304 template <class WellState>
00305 void
00306 StandardWells::
00307 computeWellConnectionDensitesPressures(const WellState& xw,
00308 const std::vector<double>& b_perf,
00309 const std::vector<double>& rsmax_perf,
00310 const std::vector<double>& rvmax_perf,
00311 const std::vector<double>& surf_dens_perf,
00312 const std::vector<double>& depth_perf,
00313 const double grav)
00314 {
00315
00316 std::vector<double> cd =
00317 WellDensitySegmented::computeConnectionDensities(
00318 wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
00319 b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
00320
00321 const int nperf = wells().well_connpos[wells().number_of_wells];
00322
00323
00324 std::vector<double> cdp =
00325 WellDensitySegmented::computeConnectionPressureDelta(
00326 wells(), depth_perf, cd, grav);
00327
00328
00329 well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf);
00330 well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf);
00331 }
00332
00333
00334
00335
00336
00337 template <class SolutionState, class WellState>
00338 void
00339 StandardWells::
00340 computeWellConnectionPressures(const SolutionState& state,
00341 const WellState& xw)
00342 {
00343 if( ! localWellsActive() ) return ;
00344
00345
00346
00347 std::vector<double> b_perf;
00348 std::vector<double> rsmax_perf;
00349 std::vector<double> rvmax_perf;
00350 std::vector<double> surf_dens_perf;
00351 computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
00352
00353 const Vector& pdepth = perf_cell_depth_;
00354 const int nperf = wells().well_connpos[wells().number_of_wells];
00355 const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
00356
00357 computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
00358
00359 }
00360
00361
00362
00363
00364
00365 template <class ReservoirResidualQuant, class SolutionState>
00366 void
00367 StandardWells::
00368 extractWellPerfProperties(const SolutionState& ,
00369 const std::vector<ReservoirResidualQuant>& rq,
00370 std::vector<ADB>& mob_perfcells,
00371 std::vector<ADB>& b_perfcells) const
00372 {
00373
00374
00375 if ( !localWellsActive() ) {
00376 mob_perfcells.clear();
00377 b_perfcells.clear();
00378 return;
00379 } else {
00380 const std::vector<int>& well_cells = wellOps().well_cells;
00381 const int num_phases = wells().number_of_phases;
00382 mob_perfcells.resize(num_phases, ADB::null());
00383 b_perfcells.resize(num_phases, ADB::null());
00384 for (int phase = 0; phase < num_phases; ++phase) {
00385 mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
00386 b_perfcells[phase] = subset(rq[phase].b, well_cells);
00387 }
00388 }
00389 }
00390
00391
00392
00393
00394
00395
00396 template <class SolutionState>
00397 void
00398 StandardWells::
00399 computeWellFlux(const SolutionState& state,
00400 const std::vector<ADB>& mob_perfcells,
00401 const std::vector<ADB>& b_perfcells,
00402 Vector& aliveWells,
00403 std::vector<ADB>& cq_s) const
00404 {
00405 if( ! localWellsActive() ) return ;
00406
00407 const int np = wells().number_of_phases;
00408 const int nw = wells().number_of_wells;
00409 const int nperf = wells().well_connpos[nw];
00410 Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
00411 const std::vector<int>& well_cells = wellOps().well_cells;
00412
00413
00414 const Vector& cdp = wellPerforationPressureDiffs();
00415
00416 const ADB& p_perfcells = subset(state.pressure, well_cells);
00417
00418
00419 const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
00420
00421
00422 const ADB drawdown = p_perfcells - perfpressure;
00423
00424
00425
00426
00427
00428 Vector selectInjectingPerforations = Vector::Zero(nperf);
00429
00430 Vector selectProducingPerforations = Vector::Zero(nperf);
00431 for (int c = 0; c < nperf; ++c){
00432 if (drawdown.value()[c] < 0)
00433 selectInjectingPerforations[c] = 1;
00434 else
00435 selectProducingPerforations[c] = 1;
00436 }
00437
00438
00439 const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
00440 const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
00441 for (int w = 0; w < nw; ++w) {
00442 if (!wells().allow_cf[w]) {
00443 for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
00444
00445
00446
00447
00448 if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
00449 selectProducingPerforations[perf] = 0.0;
00450 } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
00451 selectInjectingPerforations[perf] = 0.0;
00452 }
00453 }
00454 }
00455 }
00456
00457
00458
00459 std::vector<ADB> cq_p(np, ADB::null());
00460 std::vector<ADB> cq_ps(np, ADB::null());
00461 for (int phase = 0; phase < np; ++phase) {
00462 cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
00463 cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
00464 }
00465 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
00466 if ((*active_)[Oil] && (*active_)[Gas]) {
00467 const int oilpos = pu.phase_pos[Oil];
00468 const int gaspos = pu.phase_pos[Gas];
00469 const ADB cq_psOil = cq_ps[oilpos];
00470 const ADB cq_psGas = cq_ps[gaspos];
00471 const ADB& rv_perfcells = subset(state.rv, well_cells);
00472 const ADB& rs_perfcells = subset(state.rs, well_cells);
00473 cq_ps[gaspos] += rs_perfcells * cq_psOil;
00474 cq_ps[oilpos] += rv_perfcells * cq_psGas;
00475 }
00476
00477
00478
00479 ADB total_mob = mob_perfcells[0];
00480 for (int phase = 1; phase < np; ++phase) {
00481 total_mob += mob_perfcells[phase];
00482 }
00483
00484 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
00485
00486
00487 if (store_well_perforation_fluxes_) {
00488
00489 Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
00490 wf = cqt_i.value();
00491 for (int phase = 0; phase < np; ++phase) {
00492 wf += cq_p[phase].value();
00493 }
00494 }
00495
00496
00497
00498
00499
00500
00501 const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
00502 std::vector<ADB> wbq(np, ADB::null());
00503 ADB wbqt = ADB::constant(Vector::Zero(nw));
00504 for (int phase = 0; phase < np; ++phase) {
00505 const ADB& q_ps = wellOps().p2w * cq_ps[phase];
00506 const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
00507 Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
00508 const int pos = pu.phase_pos[phase];
00509 wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
00510 wbqt += wbq[phase];
00511 }
00512
00513 Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
00514 std::vector<ADB> cmix_s(np, ADB::null());
00515 for (int phase = 0; phase < np; ++phase) {
00516 const int pos = pu.phase_pos[phase];
00517 cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
00518 }
00519
00520
00521 ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
00522
00523 if ((*active_)[Water]) {
00524 const int watpos = pu.phase_pos[Water];
00525 volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
00526 }
00527
00528 if ((*active_)[Oil] && (*active_)[Gas]) {
00529
00530 const ADB& rv_perfcells = subset(state.rv, well_cells);
00531 const ADB& rs_perfcells = subset(state.rs, well_cells);
00532 const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
00533
00534 const int oilpos = pu.phase_pos[Oil];
00535 const int gaspos = pu.phase_pos[Gas];
00536
00537 const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * cmix_s[gaspos]) / d;
00538 volumeRatio += tmp_oil / b_perfcells[oilpos];
00539
00540 const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
00541 volumeRatio += tmp_gas / b_perfcells[gaspos];
00542 }
00543 else {
00544 if ((*active_)[Oil]) {
00545 const int oilpos = pu.phase_pos[Oil];
00546 volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
00547 }
00548 if ((*active_)[Gas]) {
00549 const int gaspos = pu.phase_pos[Gas];
00550 volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
00551 }
00552 }
00553
00554
00555
00556 ADB cqt_is = cqt_i/volumeRatio;
00557
00558
00559 cq_s.resize(np, ADB::null());
00560 for (int phase = 0; phase < np; ++phase) {
00561 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
00562 }
00563
00564
00565 aliveWells = Vector::Constant(nw, 1.0);
00566 for (int w = 0; w < nw; ++w) {
00567 if (wbqt.value()[w] == 0) {
00568 aliveWells[w] = 0.0;
00569 }
00570 }
00571 }
00572
00573
00574
00575
00576
00577 template <class SolutionState, class WellState>
00578 void
00579 StandardWells::
00580 updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
00581 const SolutionState& state,
00582 WellState& xw) const
00583 {
00584 if ( !localWellsActive() )
00585 {
00586
00587
00588 return;
00589 }
00590
00591
00592 const int np = wells().number_of_phases;
00593 const int nw = wells().number_of_wells;
00594 const int nperf = wells().well_connpos[nw];
00595 Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
00596 for (int phase = 1; phase < np; ++phase) {
00597 cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
00598 }
00599 xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
00600
00601
00602 const Vector& cdp = wellPerforationPressureDiffs();
00603 const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp;
00604 xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
00605 }
00606
00607
00608
00609
00610
00611 template <class WellState>
00612 void
00613 StandardWells::
00614 updateWellState(const Vector& dwells,
00615 const double dpmaxrel,
00616 WellState& well_state)
00617 {
00618 if( localWellsActive() )
00619 {
00620 const int np = wells().number_of_phases;
00621 const int nw = wells().number_of_wells;
00622
00623
00624 int varstart = 0;
00625 const Vector dqs = subset(dwells, Span(np*nw, 1, varstart));
00626 varstart += dqs.size();
00627 const Vector dbhp = subset(dwells, Span(nw, 1, varstart));
00628 varstart += dbhp.size();
00629 assert(varstart == dwells.size());
00630
00631
00632
00633
00634
00635
00636 const DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
00637 const Vector dwr = Eigen::Map<const Vector>(wwr.data(), nw*np);
00638 const Vector wr_old = Eigen::Map<const Vector>(&well_state.wellRates()[0], nw*np);
00639 const Vector wr = wr_old - dwr;
00640 std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
00641
00642
00643 const Vector bhp_old = Eigen::Map<const Vector>(&well_state.bhp()[0], nw, 1);
00644 const Vector dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
00645 const Vector bhp = bhp_old - dbhp_limited;
00646 std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
00647
00648
00649 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
00650
00651 #pragma omp parallel for schedule(static)
00652 for (int w = 0; w < nw; ++w) {
00653 const WellControls* wc = wells().ctrls[w];
00654 const int nwc = well_controls_get_num(wc);
00655
00656
00657
00658 for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
00659 if (well_controls_iget_type(wc, ctrl_index) == THP) {
00660 double aqua = 0.0;
00661 double liquid = 0.0;
00662 double vapour = 0.0;
00663
00664 if ((*active_)[ Water ]) {
00665 aqua = wr[w*np + pu.phase_pos[ Water ] ];
00666 }
00667 if ((*active_)[ Oil ]) {
00668 liquid = wr[w*np + pu.phase_pos[ Oil ] ];
00669 }
00670 if ((*active_)[ Gas ]) {
00671 vapour = wr[w*np + pu.phase_pos[ Gas ] ];
00672 }
00673
00674 double alq = well_controls_iget_alq(wc, ctrl_index);
00675 int table_id = well_controls_iget_vfp(wc, ctrl_index);
00676
00677 const WellType& well_type = wells().type[w];
00678 const int perf = wells().well_connpos[w];
00679 if (well_type == INJECTOR) {
00680 double dp = wellhelpers::computeHydrostaticCorrection(
00681 wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
00682 wellPerforationDensities()[perf], gravity_);
00683
00684 well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
00685 }
00686 else if (well_type == PRODUCER) {
00687 double dp = wellhelpers::computeHydrostaticCorrection(
00688 wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
00689 wellPerforationDensities()[perf], gravity_);
00690
00691 well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
00692 }
00693 else {
00694 OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
00695 }
00696
00697
00698 break;
00699 }
00700 }
00701 }
00702 }
00703 }
00704
00705
00706
00707
00708
00709 template <class WellState>
00710 void
00711 StandardWells::
00712 updateWellControls(WellState& xw) const
00713 {
00714 wellhelpers::WellSwitchingLogger logger;
00715
00716 if( !localWellsActive() ) return ;
00717
00718
00719
00720 const int np = wells().number_of_phases;
00721 const int nw = wells().number_of_wells;
00722 #pragma omp parallel for schedule(dynamic)
00723 for (int w = 0; w < nw; ++w) {
00724 const WellControls* wc = wells().ctrls[w];
00725
00726
00727
00728 int current = xw.currentControls()[w];
00729
00730
00731
00732
00733 const int nwc = well_controls_get_num(wc);
00734
00735
00736 assert(nwc != 0);
00737
00738 bool constraint_violated = false;
00739 int number_iterations = 0;
00740 const int max_iterations = 2 * nwc;
00741 do {
00742 updateWellStateWithTarget(wc, current, w, xw);
00743 int ctrl_index = 0;
00744 for (; ctrl_index < nwc; ++ctrl_index) {
00745 if (ctrl_index == current) {
00746
00747
00748
00749 continue;
00750 }
00751 if (wellhelpers::constraintBroken(
00752 xw.bhp(), xw.thp(), xw.wellRates(),
00753 w, np, wells().type[w], wc, ctrl_index)) {
00754
00755 break;
00756 }
00757 }
00758
00759
00760 if (ctrl_index != nwc) {
00761
00762
00763
00764 logger.wellSwitched(wells().name[w],
00765 well_controls_iget_type(wc, current),
00766 well_controls_iget_type(wc, ctrl_index));
00767
00768 xw.currentControls()[w] = ctrl_index;
00769 current = xw.currentControls()[w];
00770 constraint_violated = true;
00771 } else {
00772 constraint_violated = false;
00773 }
00774 ++number_iterations;
00775
00776 if (number_iterations > max_iterations) {
00777 OPM_THROW(Opm::NumericalProblem, "Could not find proper control within " << number_iterations << " iterations!");
00778 break;
00779 }
00780 } while (constraint_violated);
00781
00782
00783 if (wellCollection()->groupControlActive()) {
00784
00785
00786 WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
00787
00788
00789 if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
00790
00791 well_node.setIndividualControl(false);
00792 } else {
00793
00794 well_node.setIndividualControl(true);
00795 }
00796 }
00797 }
00798
00799 }
00800
00801
00802
00803
00804
00805 template <class SolutionState>
00806 void
00807 StandardWells::
00808 addWellFluxEq(const std::vector<ADB>& cq_s,
00809 const SolutionState& state,
00810 LinearisedBlackoilResidual& residual)
00811 {
00812 if( !localWellsActive() )
00813 {
00814
00815
00816 return;
00817 }
00818
00819 const int np = wells().number_of_phases;
00820 const int nw = wells().number_of_wells;
00821 ADB qs = state.qs;
00822 for (int phase = 0; phase < np; ++phase) {
00823 qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
00824
00825 }
00826
00827 residual.well_flux_eq = qs;
00828 }
00829
00830
00831
00832
00833
00834 template <class SolutionState, class WellState>
00835 void
00836 StandardWells::addWellControlEq(const SolutionState& state,
00837 const WellState& xw,
00838 const Vector& aliveWells,
00839 LinearisedBlackoilResidual& residual)
00840 {
00841 if( ! localWellsActive() ) return;
00842
00843 const int np = wells().number_of_phases;
00844 const int nw = wells().number_of_wells;
00845
00846 ADB aqua = ADB::constant(Vector::Zero(nw));
00847 ADB liquid = ADB::constant(Vector::Zero(nw));
00848 ADB vapour = ADB::constant(Vector::Zero(nw));
00849
00850 if ((*active_)[Water]) {
00851 aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
00852 }
00853 if ((*active_)[Oil]) {
00854 liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
00855 }
00856 if ((*active_)[Gas]) {
00857 vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
00858 }
00859
00860
00861 std::vector<int> inj_table_id(nw, -1);
00862 std::vector<int> prod_table_id(nw, -1);
00863 Vector thp_inj_target_v = Vector::Zero(nw);
00864 Vector thp_prod_target_v = Vector::Zero(nw);
00865 Vector alq_v = Vector::Zero(nw);
00866
00867
00868 Vector rho_v = Vector::Zero(nw);
00869 Vector vfp_ref_depth_v = Vector::Zero(nw);
00870
00871
00872 Vector bhp_targets = Vector::Zero(nw);
00873 Vector rate_targets = Vector::Zero(nw);
00874 Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
00875
00876
00877 std::vector<int> bhp_elems;
00878 std::vector<int> thp_inj_elems;
00879 std::vector<int> thp_prod_elems;
00880 std::vector<int> rate_elems;
00881
00882
00883
00884 for (int w = 0; w < nw; ++w) {
00885 auto wc = wells().ctrls[w];
00886
00887
00888
00889
00890 const int current = xw.currentControls()[w];
00891
00892 switch (well_controls_iget_type(wc, current)) {
00893 case BHP:
00894 {
00895 bhp_elems.push_back(w);
00896 bhp_targets(w) = well_controls_iget_target(wc, current);
00897 rate_targets(w) = -1e100;
00898 }
00899 break;
00900
00901 case THP:
00902 {
00903 const int perf = wells().well_connpos[w];
00904 rho_v[w] = wellPerforationDensities()[perf];
00905
00906 const int table_id = well_controls_iget_vfp(wc, current);
00907 const double target = well_controls_iget_target(wc, current);
00908
00909 const WellType& well_type = wells().type[w];
00910 if (well_type == INJECTOR) {
00911 inj_table_id[w] = table_id;
00912 thp_inj_target_v[w] = target;
00913 alq_v[w] = -1e100;
00914
00915 vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
00916
00917 thp_inj_elems.push_back(w);
00918 }
00919 else if (well_type == PRODUCER) {
00920 prod_table_id[w] = table_id;
00921 thp_prod_target_v[w] = target;
00922 alq_v[w] = well_controls_iget_alq(wc, current);
00923
00924 vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
00925
00926 thp_prod_elems.push_back(w);
00927 }
00928 else {
00929 OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well");
00930 }
00931 bhp_targets(w) = -1e100;
00932 rate_targets(w) = -1e100;
00933 }
00934 break;
00935
00936 case RESERVOIR_RATE:
00937 case SURFACE_RATE:
00938 {
00939 rate_elems.push_back(w);
00940
00941
00942
00943
00944 const double* const distr =
00945 well_controls_iget_distr(wc, current);
00946
00947 for (int p = 0; p < np; ++p) {
00948 rate_distr.insert(w, p*nw + w) = distr[p];
00949 }
00950
00951 bhp_targets(w) = -1.0e100;
00952 rate_targets(w) = well_controls_iget_target(wc, current);
00953 }
00954 break;
00955 }
00956 }
00957
00958
00959 const ADB thp_inj_target = ADB::constant(thp_inj_target_v);
00960 const ADB thp_prod_target = ADB::constant(thp_prod_target_v);
00961 const ADB alq = ADB::constant(alq_v);
00962 const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
00963 const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
00964
00965
00966 const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity_);
00967 const ADB dp = ADB::constant(dp_v);
00968 const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
00969 const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
00970
00971
00972 const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj;
00973 const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod;
00974 const ADB bhp_residual = state.bhp - bhp_targets;
00975 const ADB rate_residual = rate_distr * state.qs - rate_targets;
00976
00977
00978 residual.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) +
00979 superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) +
00980 superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) +
00981 superset(subset(rate_residual, rate_elems), rate_elems, nw);
00982
00983
00984
00985
00986 Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
00987 for (int w = 0; w < nw; ++w) {
00988 for (int phase = 0; phase < np; ++phase) {
00989 rate_summer.insert(w, phase*nw + w) = 1.0;
00990 }
00991 }
00992 Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
00993 residual.well_eq = alive_selector.select(residual.well_eq, rate_summer * state.qs);
00994
00995 }
00996
00997
00998
00999
01000
01001 template <class SolutionState, class WellState>
01002 void
01003 StandardWells::computeWellPotentials(const std::vector<ADB>& mob_perfcells,
01004 const std::vector<ADB>& b_perfcells,
01005 const WellState& well_state,
01006 SolutionState& state0,
01007 std::vector<double>& well_potentials) const
01008 {
01009 const int nw = wells().number_of_wells;
01010 const int np = wells().number_of_phases;
01011 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
01012
01013 Vector bhps = Vector::Zero(nw);
01014 for (int w = 0; w < nw; ++w) {
01015 const WellControls* ctrl = wells().ctrls[w];
01016 const int nwc = well_controls_get_num(ctrl);
01017
01018
01019
01020 for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
01021
01022 if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
01023 bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
01024 }
01025
01026 if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
01027 double aqua = 0.0;
01028 double liquid = 0.0;
01029 double vapour = 0.0;
01030
01031 if ((*active_)[ Water ]) {
01032 aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
01033 }
01034 if ((*active_)[ Oil ]) {
01035 liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
01036 }
01037 if ((*active_)[ Gas ]) {
01038 vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
01039 }
01040
01041 const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
01042 const double& thp = well_controls_iget_target(ctrl, ctrl_index);
01043 const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
01044
01045
01046 const WellType& well_type = wells().type[w];
01047 const int perf = wells().well_connpos[w];
01048
01049 if (well_type == INJECTOR) {
01050 double dp = wellhelpers::computeHydrostaticCorrection(
01051 wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
01052 wellPerforationDensities()[perf], gravity_);
01053 const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
01054
01055 if ( bhp < bhps[w]) {
01056 bhps[w] = bhp;
01057 }
01058 }
01059 else if (well_type == PRODUCER) {
01060 double dp = wellhelpers::computeHydrostaticCorrection(
01061 wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
01062 wellPerforationDensities()[perf], gravity_);
01063
01064 const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
01065
01066 if ( bhp > bhps[w]) {
01067 bhps[w] = bhp;
01068 }
01069 }
01070 else {
01071 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
01072 }
01073 }
01074 }
01075
01076 }
01077
01078
01079 state0.bhp = ADB::constant(bhps);
01080
01081
01082 Vector aliveWells;
01083 std::vector<ADB> perf_potentials;
01084 computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, perf_potentials);
01085
01086 well_potentials.resize(nw * np, 0.0);
01087
01088 for (int p = 0; p < np; ++p) {
01089 for (int w = 0; w < nw; ++w) {
01090 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) {
01091 well_potentials[w * np + p] += perf_potentials[p].value()[perf];
01092 }
01093 }
01094 }
01095 }
01096
01097
01098
01099
01100
01101 void
01102 StandardWells::variableStateWellIndices(std::vector<int>& indices,
01103 int& next) const
01104 {
01105 indices[Qs] = next++;
01106 indices[Bhp] = next++;
01107 }
01108
01109
01110
01111
01112
01113 template <class SolutionState>
01114 void
01115 StandardWells::
01116 variableStateExtractWellsVars(const std::vector<int>& indices,
01117 std::vector<ADB>& vars,
01118 SolutionState& state) const
01119 {
01120
01121 state.qs = std::move(vars[indices[Qs]]);
01122
01123
01124 state.bhp = std::move(vars[indices[Bhp]]);
01125 }
01126
01127
01128
01129
01130
01131 std::vector<int>
01132 StandardWells::variableWellStateIndices() const
01133 {
01134
01135
01136 std::vector<int> indices(5, -1);
01137 int next = 0;
01138
01139 variableStateWellIndices(indices, next);
01140
01141 assert(next == 2);
01142 return indices;
01143 }
01144
01145
01146
01147
01148
01149 template <class WellState>
01150 void
01151 StandardWells::variableWellStateInitials(const WellState& xw,
01152 std::vector<Vector>& vars0) const
01153 {
01154
01155 if ( localWellsActive() )
01156 {
01157
01158
01159 const int nw = wells().number_of_wells;
01160 const int np = wells().number_of_phases;
01161
01162
01163 const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
01164 const Vector qs = Eigen::Map<const V>(wrates.data(), nw*np);
01165 vars0.push_back(qs);
01166
01167
01168 assert (not xw.bhp().empty());
01169 const Vector bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
01170 vars0.push_back(bhp);
01171 }
01172 else
01173 {
01174
01175 vars0.push_back(V());
01176 vars0.push_back(V());
01177 }
01178 }
01179
01180
01181
01182
01183
01184 void
01185 StandardWells::setStoreWellPerforationFluxesFlag(const bool store_fluxes)
01186 {
01187 store_well_perforation_fluxes_ = store_fluxes;
01188 }
01189
01190
01191
01192
01193
01194 const StandardWells::Vector&
01195 StandardWells::getStoredWellPerforationFluxes() const
01196 {
01197 assert(store_well_perforation_fluxes_);
01198 return well_perforation_fluxes_;
01199 }
01200
01201
01202
01203
01204
01205
01206 template<class WellState>
01207 void
01208 StandardWells::
01209 updateListEconLimited(const Schedule& schedule,
01210 const int current_step,
01211 const Wells* wells_struct,
01212 const WellState& well_state,
01213 DynamicListEconLimited& list_econ_limited) const
01214 {
01215
01216 const int nw = ( wells_struct ) ? wells_struct->number_of_wells : 0;
01217
01218 for (int w = 0; w < nw; ++w) {
01219
01220 bool rate_limit_violated = false;
01221 const std::string& well_name = wells_struct->name[w];
01222 const Well* well_ecl = schedule.getWell(well_name);
01223 const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step);
01224
01225
01226 if (wells_struct->type[w] != PRODUCER) {
01227 continue;
01228 }
01229
01230
01231 if ( !econ_production_limits.onAnyEffectiveLimit() ) {
01232 continue;
01233 }
01234
01235
01236 const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
01237 if (quantity_limit == WellEcon::POTN) {
01238 const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n")
01239 + std::string("All the limits will be evaluated based on RATE. ");
01240 OpmLog::warning("NOT_SUPPORTING_POTN", msg);
01241 }
01242
01243 const WellMapType& well_map = well_state.wellMap();
01244 const typename WellMapType::const_iterator i_well = well_map.find(well_name);
01245 assert(i_well != well_map.end());
01246 const WellMapEntryType& map_entry = i_well->second;
01247 const int well_number = map_entry[0];
01248
01249 if (econ_production_limits.onAnyRateLimit()) {
01250 rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number);
01251 }
01252
01253 if (rate_limit_violated) {
01254 if (econ_production_limits.endRun()) {
01255 const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n")
01256 + std::string("the program will keep running after ") + well_name + std::string(" is closed");
01257 OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message);
01258 }
01259
01260 if (econ_production_limits.validFollowonWell()) {
01261 OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
01262 }
01263
01264 if (well_ecl->getAutomaticShutIn()) {
01265 list_econ_limited.addShutWell(well_name);
01266 const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit");
01267 OpmLog::info(msg);
01268 } else {
01269 list_econ_limited.addStoppedWell(well_name);
01270 const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit");
01271 OpmLog::info(msg);
01272 }
01273
01274 continue;
01275 }
01276
01277
01278 bool ratio_limits_violated = false;
01279 RatioCheckTuple ratio_check_return;
01280
01281 if (econ_production_limits.onAnyRatioLimit()) {
01282 ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry);
01283 ratio_limits_violated = std::get<0>(ratio_check_return);
01284 }
01285
01286 if (ratio_limits_violated) {
01287 const bool last_connection = std::get<1>(ratio_check_return);
01288 const int worst_offending_connection = std::get<2>(ratio_check_return);
01289
01290 const int perf_start = map_entry[1];
01291
01292 assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2]));
01293
01294 const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection];
01295 list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
01296 const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ")
01297 + well_name + std::string(" will be closed due to economic limit");
01298 OpmLog::info(msg);
01299
01300 if (last_connection) {
01301 list_econ_limited.addShutWell(well_name);
01302 const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed");
01303 OpmLog::info(msg2);
01304 }
01305 }
01306
01307 }
01308 }
01309
01310
01311
01312
01313
01314 template <class WellState>
01315 bool
01316 StandardWells::
01317 checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
01318 const WellState& well_state,
01319 const int well_number) const
01320 {
01321 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
01322 const int np = well_state.numPhases();
01323
01324 if (econ_production_limits.onMinOilRate()) {
01325 assert((*active_)[Oil]);
01326 const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
01327 const double min_oil_rate = econ_production_limits.minOilRate();
01328 if (std::abs(oil_rate) < min_oil_rate) {
01329 return true;
01330 }
01331 }
01332
01333 if (econ_production_limits.onMinGasRate() ) {
01334 assert((*active_)[Gas]);
01335 const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ];
01336 const double min_gas_rate = econ_production_limits.minGasRate();
01337 if (std::abs(gas_rate) < min_gas_rate) {
01338 return true;
01339 }
01340 }
01341
01342 if (econ_production_limits.onMinLiquidRate() ) {
01343 assert((*active_)[Oil]);
01344 assert((*active_)[Water]);
01345 const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
01346 const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
01347 const double liquid_rate = oil_rate + water_rate;
01348 const double min_liquid_rate = econ_production_limits.minLiquidRate();
01349 if (std::abs(liquid_rate) < min_liquid_rate) {
01350 return true;
01351 }
01352 }
01353
01354 if (econ_production_limits.onMinReservoirFluidRate()) {
01355 OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
01356 }
01357
01358 return false;
01359 }
01360
01361
01362
01363
01364
01365 template <class WellState>
01366 StandardWells::RatioCheckTuple
01367 StandardWells::
01368 checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
01369 const WellState& well_state,
01370 const WellMapEntryType& map_entry) const
01371 {
01372
01373
01374
01375
01376
01377
01378
01379
01380 bool any_limit_violated = false;
01381 bool last_connection = false;
01382 int worst_offending_connection = INVALIDCONNECTION;
01383 double violation_extent = -1.0;
01384
01385 if (econ_production_limits.onMaxWaterCut()) {
01386 const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry);
01387 bool water_cut_violated = std::get<0>(water_cut_return);
01388 if (water_cut_violated) {
01389 any_limit_violated = true;
01390 const double violation_extent_water_cut = std::get<3>(water_cut_return);
01391 if (violation_extent_water_cut > violation_extent) {
01392 violation_extent = violation_extent_water_cut;
01393 worst_offending_connection = std::get<2>(water_cut_return);
01394 last_connection = std::get<1>(water_cut_return);
01395 }
01396 }
01397 }
01398
01399 if (econ_production_limits.onMaxGasOilRatio()) {
01400 OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
01401 }
01402
01403 if (econ_production_limits.onMaxWaterGasRatio()) {
01404 OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
01405 }
01406
01407 if (econ_production_limits.onMaxGasLiquidRatio()) {
01408 OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
01409 }
01410
01411 if (any_limit_violated) {
01412 assert(worst_offending_connection >=0);
01413 assert(violation_extent > 1.);
01414 }
01415
01416 return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
01417 }
01418
01419
01420
01421
01422
01423 template <class WellState>
01424 StandardWells::RatioCheckTuple
01425 StandardWells::
01426 checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
01427 const WellState& well_state,
01428 const WellMapEntryType& map_entry) const
01429 {
01430 bool water_cut_limit_violated = false;
01431 int worst_offending_connection = INVALIDCONNECTION;
01432 bool last_connection = false;
01433 double violation_extent = -1.0;
01434
01435 const int np = well_state.numPhases();
01436 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
01437 const int well_number = map_entry[0];
01438
01439 assert((*active_)[Oil]);
01440 assert((*active_)[Water]);
01441
01442 const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
01443 const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
01444 const double liquid_rate = oil_rate + water_rate;
01445 double water_cut;
01446 if (std::abs(liquid_rate) != 0.) {
01447 water_cut = water_rate / liquid_rate;
01448 } else {
01449 water_cut = 0.0;
01450 }
01451
01452 const double max_water_cut_limit = econ_production_limits.maxWaterCut();
01453 if (water_cut > max_water_cut_limit) {
01454 water_cut_limit_violated = true;
01455 }
01456
01457 if (water_cut_limit_violated) {
01458
01459 const int perf_start = map_entry[1];
01460 const int perf_number = map_entry[2];
01461
01462 std::vector<double> water_cut_perf(perf_number);
01463 for (int perf = 0; perf < perf_number; ++perf) {
01464 const int i_perf = perf_start + perf;
01465 const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
01466 const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
01467 const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
01468 if (std::abs(liquid_perf_rate) != 0.) {
01469 water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
01470 } else {
01471 water_cut_perf[perf] = 0.;
01472 }
01473 }
01474
01475 last_connection = (perf_number == 1);
01476 if (last_connection) {
01477 worst_offending_connection = 0;
01478 violation_extent = water_cut_perf[0] / max_water_cut_limit;
01479 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
01480 }
01481
01482 double max_water_cut_perf = 0.;
01483 for (int perf = 0; perf < perf_number; ++perf) {
01484 if (water_cut_perf[perf] > max_water_cut_perf) {
01485 worst_offending_connection = perf;
01486 max_water_cut_perf = water_cut_perf[perf];
01487 }
01488 }
01489
01490 assert(max_water_cut_perf != 0.);
01491 assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
01492
01493 violation_extent = max_water_cut_perf / max_water_cut_limit;
01494 }
01495
01496 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
01497 }
01498
01499
01500
01501
01502
01503 WellCollection* StandardWells::wellCollection() const
01504 {
01505 return well_collection_;
01506 }
01507
01508
01509
01510
01511
01512 void StandardWells::calculateEfficiencyFactors()
01513 {
01514 if ( !localWellsActive() ) {
01515 return;
01516 }
01517
01518 const int nw = wells_->number_of_wells;
01519
01520 Vector well_efficiency_factors = Vector::Ones(nw);
01521
01522 for (int w = 0; w < nw; ++w) {
01523 const std::string well_name = wells_->name[w];
01524 const WellNode& well_node = well_collection_->findWellNode(well_name);
01525
01526 well_efficiency_factors(w) = well_node.getAccumulativeEfficiencyFactor();
01527 }
01528
01529
01530 well_perforation_efficiency_factors_ = wellOps().w2p * well_efficiency_factors.matrix();
01531 }
01532
01533
01534
01535
01536
01537 const StandardWells::Vector&
01538 StandardWells::wellPerfEfficiencyFactors() const
01539 {
01540 return well_perforation_efficiency_factors_;
01541 }
01542
01543
01544
01545
01546
01547 template <class WellState>
01548 void
01549 StandardWells::
01550 updateWellStateWithTarget(const WellControls* wc,
01551 const int current,
01552 const int well_index,
01553 WellState& xw) const
01554 {
01555 const int np = wells().number_of_phases;
01556
01557
01558
01559 const double target = well_controls_iget_target(wc, current);
01560 const double* distr = well_controls_iget_distr(wc, current);
01561 switch (well_controls_iget_type(wc, current)) {
01562 case BHP:
01563 xw.bhp()[well_index] = target;
01564 break;
01565 case THP: {
01566 double aqua = 0.0;
01567 double liquid = 0.0;
01568 double vapour = 0.0;
01569
01570 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
01571
01572 if ((*active_)[ Water ]) {
01573 aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
01574 }
01575 if ((*active_)[ Oil ]) {
01576 liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
01577 }
01578 if ((*active_)[ Gas ]) {
01579 vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
01580 }
01581
01582 const int vfp = well_controls_iget_vfp(wc, current);
01583 const double& thp = well_controls_iget_target(wc, current);
01584 const double& alq = well_controls_iget_alq(wc, current);
01585
01586
01587 const WellType& well_type = wells().type[well_index];
01588
01589 const int perf = wells().well_connpos[well_index];
01590 const double rho = well_perforation_densities_[perf];
01591
01592 if (well_type == INJECTOR) {
01593 double dp = wellhelpers::computeHydrostaticCorrection(
01594 wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
01595 rho, gravity_);
01596
01597 xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
01598 }
01599 else if (well_type == PRODUCER) {
01600 double dp = wellhelpers::computeHydrostaticCorrection(
01601 wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
01602 rho, gravity_);
01603
01604 xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
01605 }
01606 else {
01607 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
01608 }
01609 break;
01610 }
01611
01612 case RESERVOIR_RATE:
01613
01614
01615
01616
01617 break;
01618
01619 case SURFACE_RATE:
01620
01621
01622 const WellType& well_type = wells().type[well_index];
01623 if (well_type == INJECTOR) {
01624 for (int phase = 0; phase < np; ++phase) {
01625 const double& compi = wells().comp_frac[np * well_index + phase];
01626 if (compi > 0.0) {
01627 xw.wellRates()[np * well_index + phase] = target * compi;
01628 }
01629 }
01630 } else if (well_type == PRODUCER) {
01631
01632
01633
01634
01635 int numPhasesWithTargetsUnderThisControl = 0;
01636 for (int phase = 0; phase < np; ++phase) {
01637 if (distr[phase] > 0.0) {
01638 numPhasesWithTargetsUnderThisControl += 1;
01639 }
01640 }
01641 for (int phase = 0; phase < np; ++phase) {
01642 if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
01643 xw.wellRates()[np * well_index + phase] = target * distr[phase];
01644 }
01645 }
01646 } else {
01647 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
01648 }
01649
01650 break;
01651 }
01652 }
01653
01654
01655 }