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/StandardWellsSolvent.hpp>
00023
00024
00025
00026
00027 namespace Opm
00028 {
00029
00030
00031
00032
00033
00034
00035 StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection)
00036 : Base(wells_arg, well_collection)
00037 , solvent_props_(nullptr)
00038 , solvent_pos_(-1)
00039 , has_solvent_(false)
00040 {
00041 }
00042
00043
00044
00045
00046
00047 void
00048 StandardWellsSolvent::initSolvent(const SolventPropsAdFromDeck* solvent_props,
00049 const int solvent_pos,
00050 const bool has_solvent)
00051 {
00052 solvent_props_ = solvent_props;
00053 solvent_pos_ = solvent_pos;
00054 has_solvent_ = has_solvent;
00055 }
00056
00057
00058
00059
00060
00061 template<class SolutionState, class WellState>
00062 void
00063 StandardWellsSolvent::
00064 computePropertiesForWellConnectionPressures(const SolutionState& state,
00065 const WellState& xw,
00066 std::vector<double>& b_perf,
00067 std::vector<double>& rsmax_perf,
00068 std::vector<double>& rvmax_perf,
00069 std::vector<double>& surf_dens_perf)
00070 {
00071
00072
00073
00074 const int nperf = wells().well_connpos[wells().number_of_wells];
00075 const int nw = wells().number_of_wells;
00076
00077
00078 const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
00079 Vector avg_press = perf_press*0;
00080 for (int w = 0; w < nw; ++w) {
00081 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
00082 const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
00083 const double p_avg = (perf_press[perf] + p_above)/2;
00084 avg_press[perf] = p_avg;
00085 }
00086 }
00087
00088 const std::vector<int>& well_cells = wellOps().well_cells;
00089
00090
00091 const ADB perf_temp = subset(state.temperature, well_cells);
00092
00093
00094
00095
00096 const ADB avg_press_ad = ADB::constant(avg_press);
00097 std::vector<PhasePresence> perf_cond(nperf);
00098 for (int perf = 0; perf < nperf; ++perf) {
00099 perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
00100 }
00101
00102 const PhaseUsage& pu = fluid_->phaseUsage();
00103 DataBlock b(nperf, pu.num_phases);
00104
00105 const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
00106 if (pu.phase_used[BlackoilPhases::Aqua]) {
00107 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
00108 }
00109
00110 assert((*active_)[Oil]);
00111 assert((*active_)[Gas]);
00112 const ADB perf_rv = subset(state.rv, well_cells);
00113 const ADB perf_rs = subset(state.rs, well_cells);
00114 const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
00115 if (pu.phase_used[BlackoilPhases::Liquid]) {
00116 const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
00117
00118 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
00119
00120 const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
00121 rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
00122 } else {
00123 rsmax_perf.assign(0.0, nperf);
00124 }
00125 V surf_dens_copy = superset(fluid_->surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
00126 for (int phase = 1; phase < pu.num_phases; ++phase) {
00127 if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
00128 continue;
00129 }
00130 surf_dens_copy += superset(fluid_->surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
00131 }
00132
00133 if (pu.phase_used[BlackoilPhases::Vapour]) {
00134
00135
00136
00137 Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
00138 Vector rhog = fluid_->surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
00139
00140 if (has_solvent_) {
00141
00142 const Vector bs = solvent_props_->bSolvent(avg_press_ad,well_cells).value();
00143
00144
00145
00146 const int nc = state.pressure.size();
00147
00148 const ADB zero = ADB::constant(Vector::Zero(nc));
00149 const ADB& ss = state.solvent_saturation;
00150 const ADB& sg = ((*active_)[ Gas ]
00151 ? state.saturation[ pu.phase_pos[ Gas ] ]
00152 : zero);
00153
00154 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00155 Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
00156
00157 Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);
00158
00159 Vector isProducer = Vector::Zero(nperf);
00160 Vector ones = Vector::Constant(nperf,1.0);
00161 for (int w = 0; w < nw; ++w) {
00162 if(wells().type[w] == PRODUCER) {
00163 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
00164 isProducer[perf] = 1;
00165 }
00166 }
00167 }
00168
00169 F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
00170
00171 bg = bg * (ones - F_solvent);
00172 bg = bg + F_solvent * bs;
00173
00174 const Vector& rhos = solvent_props_->solventSurfaceDensity(well_cells);
00175 rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
00176 }
00177 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
00178 surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
00179
00180
00181 const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
00182 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
00183 } else {
00184 rvmax_perf.assign(0.0, nperf);
00185 }
00186
00187
00188 b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
00189 surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
00190 }
00191
00192 template <class SolutionState>
00193 void
00194 StandardWellsSolvent::
00195 computeWellFlux(const SolutionState& state,
00196 const std::vector<ADB>& mob_perfcells,
00197 const std::vector<ADB>& b_perfcells,
00198 Vector& aliveWells,
00199 std::vector<ADB>& cq_s) const
00200 {
00201 if( ! localWellsActive() ) return ;
00202
00203 const int np = wells().number_of_phases;
00204 const int nw = wells().number_of_wells;
00205 const int nperf = wells().well_connpos[nw];
00206 Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
00207 const std::vector<int>& well_cells = wellOps().well_cells;
00208
00209
00210 const Vector& cdp = wellPerforationPressureDiffs();
00211
00212 const ADB& p_perfcells = subset(state.pressure, well_cells);
00213
00214
00215 const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
00216
00217
00218 const ADB drawdown = p_perfcells - perfpressure;
00219
00220
00221
00222
00223
00224 Vector selectInjectingPerforations = Vector::Zero(nperf);
00225
00226 Vector selectProducingPerforations = Vector::Zero(nperf);
00227 for (int c = 0; c < nperf; ++c){
00228 if (drawdown.value()[c] < 0)
00229 selectInjectingPerforations[c] = 1;
00230 else
00231 selectProducingPerforations[c] = 1;
00232 }
00233
00234
00235 const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
00236 const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
00237 for (int w = 0; w < nw; ++w) {
00238 if (!wells().allow_cf[w]) {
00239 for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
00240
00241
00242
00243
00244 if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
00245 selectProducingPerforations[perf] = 0.0;
00246 } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
00247 selectInjectingPerforations[perf] = 0.0;
00248 }
00249 }
00250 }
00251 }
00252
00253
00254
00255 std::vector<ADB> cq_p(np, ADB::null());
00256 std::vector<ADB> cq_ps(np, ADB::null());
00257 for (int phase = 0; phase < np; ++phase) {
00258 cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
00259 cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
00260 }
00261
00262 Vector ones = Vector::Constant(nperf,1.0);
00263 ADB F_gas = ADB::constant(ones);
00264
00265 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
00266 if ((*active_)[Oil] && (*active_)[Gas]) {
00267 const int oilpos = pu.phase_pos[Oil];
00268 const int gaspos = pu.phase_pos[Gas];
00269 const ADB cq_psOil = cq_ps[oilpos];
00270 ADB cq_psGas = cq_ps[gaspos];
00271 const ADB& rv_perfcells = subset(state.rv, well_cells);
00272 const ADB& rs_perfcells = subset(state.rs, well_cells);
00273 cq_ps[gaspos] += rs_perfcells * cq_psOil;
00274
00275 if(has_solvent_) {
00276
00277
00278 const ADB& ss = state.solvent_saturation;
00279 const ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
00280
00281 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00282 F_gas -= subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
00283 cq_psGas = cq_psGas * F_gas;
00284 }
00285 cq_ps[oilpos] += rv_perfcells * cq_psGas;
00286 }
00287
00288
00289
00290 ADB total_mob = mob_perfcells[0];
00291 for (int phase = 1; phase < np; ++phase) {
00292 total_mob += mob_perfcells[phase];
00293 }
00294
00295 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
00296
00297
00298 if (store_well_perforation_fluxes_) {
00299
00300 Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
00301 wf = cqt_i.value();
00302 for (int phase = 0; phase < np; ++phase) {
00303 wf += cq_p[phase].value();
00304 }
00305 }
00306
00307
00308
00309
00310
00311
00312 const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
00313 std::vector<ADB> wbq(np, ADB::null());
00314 ADB wbqt = ADB::constant(Vector::Zero(nw));
00315 for (int phase = 0; phase < np; ++phase) {
00316 const ADB& q_ps = wellOps().p2w * cq_ps[phase];
00317 const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
00318 Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
00319 const int pos = pu.phase_pos[phase];
00320 wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
00321 wbqt += wbq[phase];
00322 }
00323
00324 Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
00325 std::vector<ADB> cmix_s(np, ADB::null());
00326 for (int phase = 0; phase < np; ++phase) {
00327 const int pos = pu.phase_pos[phase];
00328 cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
00329 }
00330
00331
00332 ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
00333
00334 if ((*active_)[Water]) {
00335 const int watpos = pu.phase_pos[Water];
00336 volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
00337 }
00338
00339 if ((*active_)[Oil] && (*active_)[Gas]) {
00340
00341 const ADB& rv_perfcells = subset(state.rv, well_cells);
00342 const ADB& rs_perfcells = subset(state.rs, well_cells);
00343 const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
00344
00345 const int oilpos = pu.phase_pos[Oil];
00346 const int gaspos = pu.phase_pos[Gas];
00347
00348 const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * F_gas * cmix_s[gaspos]) / d;
00349 volumeRatio += tmp_oil / b_perfcells[oilpos];
00350
00351 const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
00352 volumeRatio += tmp_gas / b_perfcells[gaspos];
00353 }
00354 else {
00355 if ((*active_)[Oil]) {
00356 const int oilpos = pu.phase_pos[Oil];
00357 volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
00358 }
00359 if ((*active_)[Gas]) {
00360 const int gaspos = pu.phase_pos[Gas];
00361 volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
00362 }
00363 }
00364
00365
00366
00367 ADB cqt_is = cqt_i/volumeRatio;
00368
00369
00370 cq_s.resize(np, ADB::null());
00371 for (int phase = 0; phase < np; ++phase) {
00372 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
00373 }
00374
00375
00376 aliveWells = Vector::Constant(nw, 1.0);
00377 for (int w = 0; w < nw; ++w) {
00378 if (wbqt.value()[w] == 0) {
00379 aliveWells[w] = 0.0;
00380 }
00381 }
00382 }
00383
00384
00385
00386
00387
00388
00389 template <class SolutionState, class WellState>
00390 void
00391 StandardWellsSolvent::
00392 computeWellConnectionPressures(const SolutionState& state,
00393 const WellState& xw)
00394 {
00395 if( ! localWellsActive() ) return ;
00396
00397
00398
00399 std::vector<double> b_perf;
00400 std::vector<double> rsmax_perf;
00401 std::vector<double> rvmax_perf;
00402 std::vector<double> surf_dens_perf;
00403 computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
00404
00405 const Vector pdepth = perf_cell_depth_;
00406 const int nperf = wells().well_connpos[wells().number_of_wells];
00407 const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
00408
00409 computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
00410
00411 }
00412
00413
00414
00415
00416
00417 template <class ReservoirResidualQuant, class SolutionState>
00418 void
00419 StandardWellsSolvent::
00420 extractWellPerfProperties(const SolutionState& state,
00421 const std::vector<ReservoirResidualQuant>& rq,
00422 std::vector<ADB>& mob_perfcells,
00423 std::vector<ADB>& b_perfcells) const
00424 {
00425 Base::extractWellPerfProperties(state, rq, mob_perfcells, b_perfcells);
00426
00427 if (has_solvent_) {
00428 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
00429 int gas_pos = pu.phase_pos[Gas];
00430 const std::vector<int>& well_cells = wellOps().well_cells;
00431 const int nperf = well_cells.size();
00432
00433
00434
00435
00436
00437 mob_perfcells[gas_pos] += subset(rq[solvent_pos_].mob, well_cells);
00438
00439
00440 const int nc = rq[solvent_pos_].mob.size();
00441
00442 const ADB zero = ADB::constant(Vector::Zero(nc));
00443 const ADB& ss = state.solvent_saturation;
00444 const ADB& sg = ((*active_)[ Gas ]
00445 ? state.saturation[ pu.phase_pos[ Gas ] ]
00446 : zero);
00447
00448 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00449 ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
00450 Vector ones = Vector::Constant(nperf,1.0);
00451
00452 b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
00453 b_perfcells[gas_pos] += (F_solvent * subset(rq[solvent_pos_].b, well_cells));
00454 }
00455 }
00456
00457
00458 }