00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 namespace Opm
00024 {
00025 template<typename TypeTag>
00026 StandardWell<TypeTag>::
00027 StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
00028 : Base(well, time_step, wells, param)
00029 , perf_densities_(number_of_perforations_)
00030 , perf_pressure_diffs_(number_of_perforations_)
00031 , primary_variables_(numWellEq, 0.0)
00032 , primary_variables_evaluation_(numWellEq)
00033 , F0_(numWellEq)
00034 {
00035 duneB_.setBuildMode( OffDiagMatWell::row_wise );
00036 duneC_.setBuildMode( OffDiagMatWell::row_wise );
00037 invDuneD_.setBuildMode( DiagMatWell::row_wise );
00038 }
00039
00040
00041
00042
00043
00044 template<typename TypeTag>
00045 void
00046 StandardWell<TypeTag>::
00047 init(const PhaseUsage* phase_usage_arg,
00048 const std::vector<bool>* active_arg,
00049 const std::vector<double>& depth_arg,
00050 const double gravity_arg,
00051 const int num_cells)
00052 {
00053 Base::init(phase_usage_arg, active_arg,
00054 depth_arg, gravity_arg, num_cells);
00055
00056 perf_depth_.resize(number_of_perforations_, 0.);
00057 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00058 const int cell_idx = well_cells_[perf];
00059 perf_depth_[perf] = depth_arg[cell_idx];
00060 }
00061
00062
00063
00064
00065
00066 invDuneD_.setSize(1, 1, 1);
00067 duneB_.setSize(1, num_cells, number_of_perforations_);
00068 duneC_.setSize(1, num_cells, number_of_perforations_);
00069
00070 for (auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) {
00071
00072 row.insert(row.index());
00073 }
00074
00075 for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
00076 for (int perf = 0 ; perf < number_of_perforations_; ++perf) {
00077 const int cell_idx = well_cells_[perf];
00078 row.insert(cell_idx);
00079 }
00080 }
00081
00082
00083 for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) {
00084 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00085 const int cell_idx = well_cells_[perf];
00086 row.insert(cell_idx);
00087 }
00088 }
00089
00090 resWell_.resize(1);
00091
00092
00093 Bx_.resize( duneB_.N() );
00094 invDrw_.resize( invDuneD_.N() );
00095 }
00096
00097
00098
00099
00100
00101 template<typename TypeTag>
00102 void StandardWell<TypeTag>::
00103 initPrimaryVariablesEvaluation() const
00104 {
00105
00106
00107
00108 for (int eqIdx = 0; eqIdx < numComponents(); ++eqIdx) {
00109 assert( (size_t)eqIdx < primary_variables_.size() );
00110
00111 primary_variables_evaluation_[eqIdx] = 0.0;
00112 primary_variables_evaluation_[eqIdx].setValue(primary_variables_[eqIdx]);
00113 primary_variables_evaluation_[eqIdx].setDerivative(numEq + eqIdx, 1.0);
00114 }
00115 }
00116
00117
00118
00119
00120
00121 template<typename TypeTag>
00122 typename StandardWell<TypeTag>::EvalWell
00123 StandardWell<TypeTag>::
00124 getBhp() const
00125 {
00126 const WellControls* wc = well_controls_;
00127 if (well_controls_get_current_type(wc) == BHP) {
00128 EvalWell bhp = 0.0;
00129 const double target_rate = well_controls_get_current_target(wc);
00130 bhp.setValue(target_rate);
00131 return bhp;
00132 } else if (well_controls_get_current_type(wc) == THP) {
00133 const int control = well_controls_get_current(wc);
00134
00135 const Opm::PhaseUsage& pu = phaseUsage();
00136 std::vector<EvalWell> rates(3, 0.0);
00137 if (active()[ Water ]) {
00138 rates[ Water ]= getQs(pu.phase_pos[ Water]);
00139 }
00140 if (active()[ Oil ]) {
00141 rates[ Oil ] = getQs(pu.phase_pos[ Oil ]);
00142 }
00143 if (active()[ Gas ]) {
00144 rates[ Gas ] = getQs(pu.phase_pos[ Gas ]);
00145 }
00146 return calculateBhpFromThp(rates, control);
00147 }
00148
00149 return primary_variables_evaluation_[XvarWell];
00150 }
00151
00152
00153
00154
00155
00156 template<typename TypeTag>
00157 typename StandardWell<TypeTag>::EvalWell
00158 StandardWell<TypeTag>::
00159 getQs(const int comp_idx) const
00160 {
00161 EvalWell qs = 0.0;
00162
00163 const WellControls* wc = well_controls_;
00164 const int np = number_of_phases_;
00165 const double target_rate = well_controls_get_current_target(wc);
00166
00167 assert(comp_idx < numComponents());
00168 const auto pu = phaseUsage();
00169
00170
00171
00172 if (well_type_ == INJECTOR) {
00173 if (has_solvent) {
00174
00175
00176 double comp_frac = 0.0;
00177 if (has_solvent && comp_idx == contiSolventEqIdx) {
00178 comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
00179 } else if (comp_idx == pu.phase_pos[ Gas ]) {
00180 comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent());
00181 } else {
00182 comp_frac = comp_frac_[comp_idx];
00183 }
00184 if (comp_frac == 0.0) {
00185 return qs;
00186 }
00187
00188 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
00189 return comp_frac * primary_variables_evaluation_[XvarWell];
00190 }
00191
00192 qs.setValue(comp_frac * target_rate);
00193 return qs;
00194 }
00195
00196 const double comp_frac = comp_frac_[comp_idx];
00197 if (comp_frac == 0.0) {
00198 return qs;
00199 }
00200
00201 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
00202 return primary_variables_evaluation_[XvarWell];
00203 }
00204 qs.setValue(target_rate);
00205 return qs;
00206 }
00207
00208
00209 if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
00210 return primary_variables_evaluation_[XvarWell] * wellVolumeFractionScaled(comp_idx);
00211 }
00212
00213 if (well_controls_get_current_type(wc) == SURFACE_RATE) {
00214
00215
00216 const double* distr = well_controls_get_current_distr(wc);
00217 int num_phases_under_rate_control = 0;
00218 for (int phase = 0; phase < np; ++phase) {
00219 if (distr[phase] > 0.0) {
00220 num_phases_under_rate_control += 1;
00221 }
00222 }
00223
00224
00225 assert(num_phases_under_rate_control > 0);
00226
00227
00228 if (num_phases_under_rate_control == 1) {
00229
00230
00231 int phase_under_control = -1;
00232 for (int phase = 0; phase < np; ++phase) {
00233 if (distr[phase] > 0.0) {
00234 phase_under_control = phase;
00235 break;
00236 }
00237 }
00238
00239 assert(phase_under_control >= 0);
00240
00241 EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(phase_under_control);
00242 if (has_solvent && phase_under_control == Gas) {
00243
00244 wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx);
00245 }
00246
00247 if (comp_idx == phase_under_control) {
00248 if (has_solvent && phase_under_control == Gas) {
00249 qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() );
00250 return qs;
00251 }
00252 qs.setValue(target_rate);
00253 return qs;
00254 }
00255
00256
00257 const double eps = 1e-6;
00258 if (wellVolumeFractionScaledPhaseUnderControl < eps) {
00259 return qs;
00260 }
00261 return (target_rate * wellVolumeFractionScaled(comp_idx) / wellVolumeFractionScaledPhaseUnderControl);
00262 }
00263
00264
00265
00266 if (num_phases_under_rate_control == 2) {
00267 EvalWell combined_volume_fraction = 0.;
00268 for (int p = 0; p < np; ++p) {
00269 if (distr[p] == 1.0) {
00270 combined_volume_fraction += wellVolumeFractionScaled(p);
00271 }
00272 }
00273 return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction);
00274 }
00275
00276
00277 if (num_phases_under_rate_control == 3) {
00278 return target_rate * wellSurfaceVolumeFraction(comp_idx);
00279 }
00280 } else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
00281
00282 return target_rate * wellVolumeFractionScaled(comp_idx);
00283 } else {
00284 OPM_THROW(std::logic_error, "Unknown control type for well " << name());
00285 }
00286
00287
00288 return qs;
00289 }
00290
00291
00292
00293
00294
00295
00296 template<typename TypeTag>
00297 typename StandardWell<TypeTag>::EvalWell
00298 StandardWell<TypeTag>::
00299 wellVolumeFractionScaled(const int compIdx) const
00300 {
00301
00302 const double scal = scalingFactor(compIdx);
00303 if (scal > 0)
00304 return wellVolumeFraction(compIdx) / scal;
00305
00306
00307 return wellVolumeFraction(compIdx);
00308 }
00309
00310
00311
00312
00313
00314 template<typename TypeTag>
00315 typename StandardWell<TypeTag>::EvalWell
00316 StandardWell<TypeTag>::
00317 wellVolumeFraction(const int compIdx) const
00318 {
00319 const auto& pu = phaseUsage();
00320 if (active()[Water] && compIdx == pu.phase_pos[Water]) {
00321 return primary_variables_evaluation_[WFrac];
00322 }
00323
00324 if (active()[Gas] && compIdx == pu.phase_pos[Gas]) {
00325 return primary_variables_evaluation_[GFrac];
00326 }
00327
00328 if (has_solvent && compIdx == contiSolventEqIdx) {
00329 return primary_variables_evaluation_[SFrac];
00330 }
00331
00332
00333 EvalWell well_fraction = 1.0;
00334 if (active()[Water]) {
00335 well_fraction -= primary_variables_evaluation_[WFrac];
00336 }
00337
00338 if (active()[Gas]) {
00339 well_fraction -= primary_variables_evaluation_[GFrac];
00340 }
00341 if (has_solvent) {
00342 well_fraction -= primary_variables_evaluation_[SFrac];
00343 }
00344 return well_fraction;
00345 }
00346
00347
00348
00349
00350
00351 template<typename TypeTag>
00352 typename StandardWell<TypeTag>::EvalWell
00353 StandardWell<TypeTag>::
00354 wellSurfaceVolumeFraction(const int compIdx) const
00355 {
00356 EvalWell sum_volume_fraction_scaled = 0.;
00357 const int numComp = numComponents();
00358 for (int idx = 0; idx < numComp; ++idx) {
00359 sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
00360 }
00361
00362 assert(sum_volume_fraction_scaled.value() != 0.);
00363
00364 return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled;
00365 }
00366
00367
00368
00369
00370
00371 template<typename TypeTag>
00372 typename StandardWell<TypeTag>::EvalWell
00373 StandardWell<TypeTag>::
00374 extendEval(const Eval& in) const
00375 {
00376 EvalWell out = 0.0;
00377 out.setValue(in.value());
00378 for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
00379 out.setDerivative(eqIdx, in.derivative(eqIdx));
00380 }
00381 return out;
00382 }
00383
00384
00385
00386
00387
00388 template<typename TypeTag>
00389 void
00390 StandardWell<TypeTag>::
00391 computePerfRate(const IntensiveQuantities& intQuants,
00392 const std::vector<EvalWell>& mob_perfcells_dense,
00393 const double Tw, const EvalWell& bhp, const double& cdp,
00394 const bool& allow_cf, std::vector<EvalWell>& cq_s) const
00395 {
00396 const Opm::PhaseUsage& pu = phaseUsage();
00397 const int np = number_of_phases_;
00398 const int numComp = numComponents();
00399 std::vector<EvalWell> cmix_s(numComp,0.0);
00400 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
00401 cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
00402 }
00403 const auto& fs = intQuants.fluidState();
00404
00405 const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
00406 const EvalWell rs = extendEval(fs.Rs());
00407 const EvalWell rv = extendEval(fs.Rv());
00408 std::vector<EvalWell> b_perfcells_dense(numComp, 0.0);
00409 for (int phase = 0; phase < np; ++phase) {
00410 const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
00411 b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
00412 }
00413 if (has_solvent) {
00414 b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
00415 }
00416
00417
00418 const EvalWell well_pressure = bhp + cdp;
00419 const EvalWell drawdown = pressure - well_pressure;
00420
00421
00422 if ( drawdown.value() > 0 ) {
00423
00424 if (!allow_cf && well_type_ == INJECTOR) {
00425 return;
00426 }
00427
00428
00429 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
00430 const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown);
00431 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
00432 }
00433
00434 if (active()[Oil] && active()[Gas]) {
00435 const int oilpos = pu.phase_pos[Oil];
00436 const int gaspos = pu.phase_pos[Gas];
00437 const EvalWell cq_sOil = cq_s[oilpos];
00438 const EvalWell cq_sGas = cq_s[gaspos];
00439 cq_s[gaspos] += rs * cq_sOil;
00440 cq_s[oilpos] += rv * cq_sGas;
00441 }
00442
00443 } else {
00444
00445 if (!allow_cf && well_type_ == PRODUCER) {
00446 return;
00447 }
00448
00449
00450 EvalWell total_mob_dense = mob_perfcells_dense[0];
00451 for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) {
00452 total_mob_dense += mob_perfcells_dense[componentIdx];
00453 }
00454
00455
00456 const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
00457
00458
00459 EvalWell volumeRatio = 0.0;
00460 if (active()[Water]) {
00461 const int watpos = pu.phase_pos[Water];
00462 volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos];
00463 }
00464
00465 if (has_solvent) {
00466 volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx];
00467 }
00468
00469 if (active()[Oil] && active()[Gas]) {
00470 const int oilpos = pu.phase_pos[Oil];
00471 const int gaspos = pu.phase_pos[Gas];
00472
00473
00474 const EvalWell d = 1.0 - rv * rs;
00475
00476 if (d.value() == 0.0) {
00477 OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << name() << " during flux calcuation"
00478 << " with rs " << rs << " and rv " << rv);
00479 }
00480
00481 const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
00482
00483 volumeRatio += tmp_oil / b_perfcells_dense[oilpos];
00484
00485 const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
00486
00487 volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
00488 }
00489 else {
00490 if (active()[Oil]) {
00491 const int oilpos = pu.phase_pos[Oil];
00492 volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
00493 }
00494 if (active()[Gas]) {
00495 const int gaspos = pu.phase_pos[Gas];
00496 volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
00497 }
00498 }
00499
00500
00501 EvalWell cqt_is = cqt_i/volumeRatio;
00502
00503 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
00504 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
00505 }
00506 }
00507 }
00508
00509
00510
00511
00512
00513 template<typename TypeTag>
00514 void
00515 StandardWell<TypeTag>::
00516 assembleWellEq(Simulator& ebosSimulator,
00517 const double dt,
00518 WellState& well_state,
00519 bool only_wells)
00520 {
00521 const int numComp = numComponents();
00522 const int np = number_of_phases_;
00523
00524
00525 if (!only_wells) {
00526 duneB_ = 0.0;
00527 duneC_ = 0.0;
00528 }
00529 invDuneD_ = 0.0;
00530 resWell_ = 0.0;
00531
00532 auto& ebosJac = ebosSimulator.model().linearizer().matrix();
00533 auto& ebosResid = ebosSimulator.model().linearizer().residual();
00534
00535
00536 const double volume = 0.002831684659200;
00537
00538 const bool allow_cf = crossFlowAllowed(ebosSimulator);
00539
00540 const EvalWell& bhp = getBhp();
00541
00542 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00543
00544 const int cell_idx = well_cells_[perf];
00545 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
00546 std::vector<EvalWell> cq_s(numComp,0.0);
00547 std::vector<EvalWell> mob(numComp, 0.0);
00548 getMobility(ebosSimulator, perf, mob);
00549 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
00550
00551 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
00552
00553 const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_;
00554
00555 if (!only_wells) {
00556
00557
00558 ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value();
00559 }
00560
00561
00562 resWell_[0][componentIdx] -= cq_s_effective.value();
00563
00564
00565 for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
00566 if (!only_wells) {
00567
00568 duneC_[0][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq);
00569 }
00570 invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq);
00571 }
00572
00573 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
00574 if (!only_wells) {
00575
00576 ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx);
00577 duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
00578 }
00579 }
00580
00581
00582 if (has_solvent && componentIdx == contiSolventEqIdx) {
00583 well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
00584 } else {
00585 well_state.perfPhaseRates()[(first_perf_ + perf) * np + componentIdx] = cq_s[componentIdx].value();
00586 }
00587 }
00588
00589 if (has_polymer) {
00590 EvalWell cq_s_poly = cq_s[Water];
00591 if (well_type_ == INJECTOR) {
00592 cq_s_poly *= wpolymer();
00593 } else {
00594 cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
00595 }
00596 if (!only_wells) {
00597
00598 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
00599 ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
00600 }
00601 ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
00602 }
00603 }
00604
00605
00606 well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
00607 }
00608
00609
00610 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
00611 EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
00612 resWell_loc += getQs(componentIdx) * well_efficiency_factor_;
00613 for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
00614 invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
00615 }
00616 resWell_[0][componentIdx] += resWell_loc.value();
00617 }
00618
00619
00620 invDuneD_[0][0].invert();
00621 }
00622
00623
00624
00625
00626
00627 template<typename TypeTag>
00628 bool
00629 StandardWell<TypeTag>::
00630 crossFlowAllowed(const Simulator& ebosSimulator) const
00631 {
00632 if (getAllowCrossFlow()) {
00633 return true;
00634 }
00635
00636
00637
00638
00639
00640 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00641 const int cell_idx = well_cells_[perf];
00642 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
00643 const auto& fs = intQuants.fluidState();
00644 EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
00645 EvalWell bhp = getBhp();
00646
00647
00648 EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
00649 EvalWell drawdown = pressure - well_pressure;
00650
00651 if (drawdown.value() < 0 && well_type_ == INJECTOR) {
00652 return false;
00653 }
00654
00655 if (drawdown.value() > 0 && well_type_ == PRODUCER) {
00656 return false;
00657 }
00658 }
00659 return true;
00660 }
00661
00662
00663
00664
00665
00666 template<typename TypeTag>
00667 void
00668 StandardWell<TypeTag>::
00669 getMobility(const Simulator& ebosSimulator,
00670 const int perf,
00671 std::vector<EvalWell>& mob) const
00672 {
00673 const int np = number_of_phases_;
00674 const int cell_idx = well_cells_[perf];
00675 assert (int(mob.size()) == numComponents());
00676 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
00677 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
00678
00679
00680
00681 const int satid = saturation_table_number_[perf] - 1;
00682 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
00683 if( satid == satid_elem ) {
00684
00685 for (int phase = 0; phase < np; ++phase) {
00686 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
00687 mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
00688 }
00689 if (has_solvent) {
00690 mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
00691 }
00692 } else {
00693
00694 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
00695 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
00696 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
00697
00698
00699 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
00700
00701
00702 for (int phase = 0; phase < np; ++phase) {
00703 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
00704 mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
00705 }
00706
00707
00708 if (has_solvent) {
00709 OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
00710 }
00711 }
00712
00713
00714 if (has_polymer) {
00715
00716 EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
00717
00718 if (well_type_ == INJECTOR) {
00719 const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
00720 mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, true) );
00721 }
00722
00723 if (PolymerModule::hasPlyshlog()) {
00724
00725 const int numComp = numComponents();
00726 const bool allow_cf = crossFlowAllowed(ebosSimulator);
00727 const EvalWell& bhp = getBhp();
00728 std::vector<EvalWell> cq_s(numComp,0.0);
00729 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
00730
00731 double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];
00732 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
00733 const auto& scaledDrainageInfo =
00734 materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
00735 const Scalar& Swcr = scaledDrainageInfo.Swcr;
00736 const EvalWell poro = extendEval(intQuants.porosity());
00737 const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
00738
00739 const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
00740 EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
00741
00742 if (PolymerModule::hasShrate()) {
00743
00744
00745
00746 waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / bore_diameters_[perf];
00747 }
00748 EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
00749 EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
00750 intQuants.pvtRegionIndex(),
00751 waterVelocity);
00752
00753
00754 mob[ Water ] /= shearFactor;
00755 }
00756 }
00757 }
00758
00759
00760
00761
00762
00763 template<typename TypeTag>
00764 void
00765 StandardWell<TypeTag>::
00766 updateWellState(const BVectorWell& dwells,
00767 WellState& well_state) const
00768 {
00769 const int np = number_of_phases_;
00770 const double dBHPLimit = param_.dbhp_max_rel_;
00771 const double dFLimit = param_.dwell_fraction_max_;
00772 const auto pu = phaseUsage();
00773
00774 const std::vector<double> xvar_well_old = primary_variables_;
00775
00776
00777 std::vector<double> F(np,0.0);
00778 if (active()[ Water ]) {
00779 const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
00780 const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit);
00781 primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited;
00782 }
00783
00784 if (active()[ Gas ]) {
00785 const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
00786 const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit);
00787 primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited;
00788 }
00789
00790 if (has_solvent) {
00791 const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
00792 const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]),dFLimit);
00793 primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited;
00794 }
00795
00796 assert(active()[ Oil ]);
00797 F[pu.phase_pos[Oil]] = 1.0;
00798
00799 if (active()[ Water ]) {
00800 F[pu.phase_pos[Water]] = primary_variables_[WFrac];
00801 F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
00802 }
00803
00804 if (active()[ Gas ]) {
00805 F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
00806 F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
00807 }
00808
00809 double F_solvent = 0.0;
00810 if (has_solvent) {
00811 F_solvent = primary_variables_[SFrac];
00812 F[pu.phase_pos[Oil]] -= F_solvent;
00813 }
00814
00815 if (active()[ Water ]) {
00816 if (F[Water] < 0.0) {
00817 if (active()[ Gas ]) {
00818 F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
00819 }
00820 if (has_solvent) {
00821 F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
00822 }
00823 F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
00824 F[pu.phase_pos[Water]] = 0.0;
00825 }
00826 }
00827
00828 if (active()[ Gas ]) {
00829 if (F[pu.phase_pos[Gas]] < 0.0) {
00830 if (active()[ Water ]) {
00831 F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
00832 }
00833 if (has_solvent) {
00834 F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
00835 }
00836 F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
00837 F[pu.phase_pos[Gas]] = 0.0;
00838 }
00839 }
00840
00841 if (F[pu.phase_pos[Oil]] < 0.0) {
00842 if (active()[ Water ]) {
00843 F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
00844 }
00845 if (active()[ Gas ]) {
00846 F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
00847 }
00848 if (has_solvent) {
00849 F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
00850 }
00851 F[pu.phase_pos[Oil]] = 0.0;
00852 }
00853
00854 if (active()[ Water ]) {
00855 primary_variables_[WFrac] = F[pu.phase_pos[Water]];
00856 }
00857 if (active()[ Gas ]) {
00858 primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
00859 }
00860 if(has_solvent) {
00861 primary_variables_[SFrac] = F_solvent;
00862 }
00863
00864
00865
00866 if (has_solvent){
00867 F[pu.phase_pos[Gas]] += F_solvent;
00868 }
00869
00870
00871 const WellControls* wc = well_controls_;
00872
00873
00874
00875 const int current = well_state.currentControls()[index_of_well_];
00876 const double target_rate = well_controls_iget_target(wc, current);
00877
00878 for (int p = 0; p < np; ++p) {
00879 const double scal = scalingFactor(p);
00880 if (scal > 0) {
00881 F[p] /= scal ;
00882 } else {
00883 F[p] = 0.;
00884 }
00885 }
00886
00887 switch (well_controls_iget_type(wc, current)) {
00888 case THP:
00889 case BHP:
00890 {
00891 primary_variables_[XvarWell] = xvar_well_old[XvarWell] - dwells[0][XvarWell];
00892
00893 switch (well_type_) {
00894 case INJECTOR:
00895 for (int p = 0; p < np; ++p) {
00896 const double comp_frac = comp_frac_[p];
00897 well_state.wellRates()[index_of_well_ * np + p] = comp_frac * primary_variables_[XvarWell];
00898 }
00899 break;
00900 case PRODUCER:
00901 for (int p = 0; p < np; ++p) {
00902 well_state.wellRates()[index_of_well_ * np + p] = primary_variables_[XvarWell] * F[p];
00903 }
00904 break;
00905 }
00906
00907 if (well_controls_iget_type(wc, current) == THP) {
00908
00909
00910 std::vector<double> rates(3, 0.0);
00911
00912 const Opm::PhaseUsage& pu = phaseUsage();
00913 if (active()[ Water ]) {
00914 rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
00915 }
00916 if (active()[ Oil ]) {
00917 rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
00918 }
00919 if (active()[ Gas ]) {
00920 rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
00921 }
00922
00923 well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current);
00924 }
00925 }
00926 break;
00927 case SURFACE_RATE:
00928 case RESERVOIR_RATE:
00929 {
00930 const int sign1 = dwells[0][XvarWell] > 0 ? 1: -1;
00931 const double dx1_limited = sign1 * std::min(std::abs(dwells[0][XvarWell]),std::abs(xvar_well_old[XvarWell])*dBHPLimit);
00932 primary_variables_[XvarWell] = std::max(xvar_well_old[XvarWell] - dx1_limited,1e5);
00933 well_state.bhp()[index_of_well_] = primary_variables_[XvarWell];
00934
00935 if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
00936 if (well_type_ == PRODUCER) {
00937
00938 const double* distr = well_controls_iget_distr(wc, current);
00939
00940 double F_target = 0.0;
00941 for (int p = 0; p < np; ++p) {
00942 F_target += distr[p] * F[p];
00943 }
00944 for (int p = 0; p < np; ++p) {
00945 well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate / F_target;
00946 }
00947 } else {
00948
00949 for (int p = 0; p < np; ++p) {
00950 well_state.wellRates()[index_of_well_ * np + p] = comp_frac_[p] * target_rate;
00951 }
00952 }
00953 } else {
00954 for (int p = 0; p < np; ++p) {
00955 well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate;
00956 }
00957 }
00958 }
00959 break;
00960 }
00961
00962
00963
00964
00965 const int nwc = well_controls_get_num(wc);
00966
00967 int ctrl_index = 0;
00968 for ( ; ctrl_index < nwc; ++ctrl_index) {
00969 if (well_controls_iget_type(wc, ctrl_index) == THP) {
00970
00971 const int current = well_state.currentControls()[index_of_well_];
00972
00973 if (current == ctrl_index) {
00974 const double thp_target = well_controls_iget_target(wc, current);
00975 well_state.thp()[index_of_well_] = thp_target;
00976 } else {
00977 const Opm::PhaseUsage& pu = phaseUsage();
00978 std::vector<double> rates(3, 0.0);
00979
00980 if (active()[ Water ]) {
00981 rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ];
00982 }
00983 if (active()[ Oil ]) {
00984 rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ];
00985 }
00986 if (active()[ Gas ]) {
00987 rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ];
00988 }
00989
00990 const double bhp = well_state.bhp()[index_of_well_];
00991
00992 well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp);
00993 }
00994
00995
00996 break;
00997 }
00998 }
00999
01000
01001 if (ctrl_index == nwc) {
01002 well_state.thp()[index_of_well_] = 0.0;
01003 }
01004 }
01005
01006
01007
01008
01009
01010 template<typename TypeTag>
01011 void
01012 StandardWell<TypeTag>::
01013 updateWellStateWithTarget(const int current,
01014 WellState& xw) const
01015 {
01016
01017 const int np = number_of_phases_;
01018 const int well_index = index_of_well_;
01019 const WellControls* wc = well_controls_;
01020
01021
01022 const double target = well_controls_iget_target(wc, current);
01023 const double* distr = well_controls_iget_distr(wc, current);
01024 switch (well_controls_iget_type(wc, current)) {
01025 case BHP:
01026 xw.bhp()[well_index] = target;
01027
01028
01029 break;
01030
01031 case THP: {
01032 xw.thp()[well_index] = target;
01033
01034 const Opm::PhaseUsage& pu = phaseUsage();
01035 std::vector<double> rates(3, 0.0);
01036 if (active()[ Water ]) {
01037 rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
01038 }
01039 if (active()[ Oil ]) {
01040 rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
01041 }
01042 if (active()[ Gas ]) {
01043 rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
01044 }
01045
01046 xw.bhp()[well_index] = calculateBhpFromThp(rates, current);
01047 break;
01048 }
01049
01050 case RESERVOIR_RATE:
01051 case SURFACE_RATE:
01052
01053 int numPhasesWithTargetsUnderThisControl = 0;
01054 for (int phase = 0; phase < np; ++phase) {
01055 if (distr[phase] > 0.0) {
01056 numPhasesWithTargetsUnderThisControl += 1;
01057 }
01058 }
01059
01060 assert(numPhasesWithTargetsUnderThisControl > 0);
01061
01062 if (well_type_ == INJECTOR) {
01063
01064
01065 assert(numPhasesWithTargetsUnderThisControl == 1);
01066
01067 for (int phase = 0; phase < np; ++phase) {
01068 if (distr[phase] > 0.) {
01069 xw.wellRates()[np*well_index + phase] = target / distr[phase];
01070 } else {
01071 xw.wellRates()[np * well_index + phase] = 0.;
01072 }
01073 }
01074 } else if (well_type_ == PRODUCER) {
01075
01076
01077
01078 double original_rates_under_phase_control = 0.0;
01079 for (int phase = 0; phase < np; ++phase) {
01080 if (distr[phase] > 0.0) {
01081 original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase];
01082 }
01083 }
01084
01085 if (original_rates_under_phase_control != 0.0 ) {
01086 double scaling_factor = target / original_rates_under_phase_control;
01087
01088 for (int phase = 0; phase < np; ++phase) {
01089 xw.wellRates()[np * well_index + phase] *= scaling_factor;
01090 }
01091 } else {
01092
01093 const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
01094 for (int phase = 0; phase < np; ++phase) {
01095 if (distr[phase] > 0.0) {
01096 xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
01097 } else {
01098
01099 xw.wellRates()[np * well_index + phase] = target_rate_divided;
01100 }
01101 }
01102 }
01103 } else {
01104 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
01105 }
01106
01107 break;
01108 }
01109
01110 updatePrimaryVariables(xw);
01111 }
01112
01113
01114
01115
01116
01117 template<typename TypeTag>
01118 void
01119 StandardWell<TypeTag>::
01120 computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
01121 const WellState& xw,
01122 std::vector<double>& b_perf,
01123 std::vector<double>& rsmax_perf,
01124 std::vector<double>& rvmax_perf,
01125 std::vector<double>& surf_dens_perf) const
01126 {
01127 const int nperf = number_of_perforations_;
01128
01129 const int numComp = numComponents();
01130 const PhaseUsage& pu = phaseUsage();
01131 b_perf.resize(nperf*numComp);
01132 surf_dens_perf.resize(nperf*numComp);
01133 const int w = index_of_well_;
01134
01135
01136 if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) {
01137 rsmax_perf.resize(nperf);
01138 rvmax_perf.resize(nperf);
01139 }
01140
01141
01142 for (int perf = 0; perf < nperf; ++perf) {
01143 const int cell_idx = well_cells_[perf];
01144 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
01145 const auto& fs = intQuants.fluidState();
01146
01147
01148
01149 const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1];
01150 const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2;
01151 const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
01152
01153 if (pu.phase_used[BlackoilPhases::Aqua]) {
01154 b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * numComp] =
01155 FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
01156 }
01157
01158 if (pu.phase_used[BlackoilPhases::Vapour]) {
01159 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp;
01160 const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
01161
01162 if (pu.phase_used[BlackoilPhases::Liquid]) {
01163 const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
01164 const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
01165 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
01166 if (oilrate > 0) {
01167 const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
01168 double rv = 0.0;
01169 if (gasrate > 0) {
01170 rv = oilrate / gasrate;
01171 }
01172 rv = std::min(rv, rvmax_perf[perf]);
01173
01174 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
01175 }
01176 else {
01177 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
01178 }
01179
01180 } else {
01181 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
01182 }
01183 }
01184
01185 if (pu.phase_used[BlackoilPhases::Liquid]) {
01186 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp;
01187 const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
01188 if (pu.phase_used[BlackoilPhases::Vapour]) {
01189 rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
01190 const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
01191 const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
01192 if (gasrate > 0) {
01193 const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
01194 double rs = 0.0;
01195 if (oilrate > 0) {
01196 rs = gasrate / oilrate;
01197 }
01198 rs = std::min(rs, rsmax_perf[perf]);
01199 b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
01200 } else {
01201 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
01202 }
01203 } else {
01204 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
01205 }
01206 }
01207
01208
01209 for (int p = 0; p < pu.num_phases; ++p) {
01210 surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex());
01211 }
01212
01213
01214 if (has_solvent) {
01215 b_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
01216 surf_dens_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventRefDensity();
01217 }
01218 }
01219 }
01220
01221
01222
01223
01224
01225 template<typename TypeTag>
01226 void
01227 StandardWell<TypeTag>::
01228 computeConnectionDensities(const std::vector<double>& perfComponentRates,
01229 const std::vector<double>& b_perf,
01230 const std::vector<double>& rsmax_perf,
01231 const std::vector<double>& rvmax_perf,
01232 const std::vector<double>& surf_dens_perf)
01233 {
01234
01235 const int np = number_of_phases_;
01236 const int nperf = number_of_perforations_;
01237 const int num_comp = numComponents();
01238 const PhaseUsage& phase_usage = phaseUsage();
01239
01240
01241
01242
01243
01244 std::vector<double> q_out_perf(nperf*num_comp);
01245
01246
01247
01248 for (int perf = nperf - 1; perf >= 0; --perf) {
01249 for (int component = 0; component < num_comp; ++component) {
01250 if (perf == nperf - 1) {
01251
01252 q_out_perf[perf*num_comp+ component] = 0.0;
01253 } else {
01254
01255 q_out_perf[perf*num_comp + component] = q_out_perf[(perf+1)*num_comp + component];
01256 }
01257
01258 q_out_perf[perf*num_comp + component] -= perfComponentRates[perf*num_comp + component];
01259 }
01260 }
01261
01262
01263
01264
01265
01266 const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
01267 const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
01268 std::vector<double> mix(num_comp,0.0);
01269 std::vector<double> x(num_comp);
01270 std::vector<double> surf_dens(num_comp);
01271 std::vector<double> dens(nperf);
01272
01273 for (int perf = 0; perf < nperf; ++perf) {
01274
01275 const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf,
01276 q_out_perf.begin() + num_comp*(perf+1), 0.0);
01277 if (tot_surf_rate != 0.0) {
01278 for (int component = 0; component < num_comp; ++component) {
01279 mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate);
01280 }
01281 } else {
01282
01283 for (int phase = 0; phase < np; ++phase) {
01284 mix[phase] = comp_frac_[phase];
01285 }
01286
01287 }
01288
01289 x = mix;
01290 double rs = 0.0;
01291 double rv = 0.0;
01292 if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
01293 rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
01294 }
01295 if (!rvmax_perf.empty() && mix[gaspos] > 0.0) {
01296 rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
01297 }
01298 if (rs != 0.0) {
01299
01300 x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
01301 }
01302 if (rv != 0.0) {
01303
01304 x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
01305 }
01306 double volrat = 0.0;
01307 for (int component = 0; component < num_comp; ++component) {
01308 volrat += x[component] / b_perf[perf*num_comp+ component];
01309 }
01310 for (int component = 0; component < num_comp; ++component) {
01311 surf_dens[component] = surf_dens_perf[perf*num_comp+ component];
01312 }
01313
01314
01315 perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
01316 }
01317 }
01318
01319
01320
01321
01322 template<typename TypeTag>
01323 void
01324 StandardWell<TypeTag>::
01325 computeConnectionPressureDelta()
01326 {
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341 const int nperf = number_of_perforations_;
01342 perf_pressure_diffs_.resize(nperf, 0.0);
01343
01344 for (int perf = 0; perf < nperf; ++perf) {
01345 const double z_above = perf == 0 ? ref_depth_ : perf_depth_[perf - 1];
01346 const double dz = perf_depth_[perf] - z_above;
01347 perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * gravity_;
01348 }
01349
01350
01351
01352
01353
01354 const auto beg = perf_pressure_diffs_.begin();
01355 const auto end = perf_pressure_diffs_.end();
01356 std::partial_sum(beg, end, beg);
01357 }
01358
01359
01360
01361
01362
01363 template<typename TypeTag>
01364 typename StandardWell<TypeTag>::ConvergenceReport
01365 StandardWell<TypeTag>::
01366 getWellConvergence(const std::vector<double>& B_avg) const
01367 {
01368 typedef double Scalar;
01369 typedef std::vector< Scalar > Vector;
01370
01371 const int np = number_of_phases_;
01372 const int numComp = numComponents();
01373
01374
01375
01376 assert((int(B_avg.size()) == numComp) || has_polymer);
01377
01378 const double tol_wells = param_.tolerance_wells_;
01379 const double maxResidualAllowed = param_.max_residual_allowed_;
01380
01381
01382
01383 std::vector<Scalar> res(numComp);
01384 for (int comp = 0; comp < numComp; ++comp) {
01385
01386 res[comp] = std::abs(resWell_[0][comp]);
01387 }
01388
01389 Vector well_flux_residual(numComp);
01390
01391
01392 for ( int compIdx = 0; compIdx < numComp; ++compIdx )
01393 {
01394 well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx];
01395 }
01396
01397 ConvergenceReport report;
01398
01399
01400 for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
01401 const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
01402
01403 if (std::isnan(well_flux_residual[phaseIdx])) {
01404 report.nan_residual_found = true;
01405 const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
01406 report.nan_residual_wells.push_back(problem_well);
01407 } else {
01408 if (well_flux_residual[phaseIdx] > maxResidualAllowed) {
01409 report.too_large_residual_found = true;
01410 const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
01411 report.too_large_residual_wells.push_back(problem_well);
01412 }
01413 }
01414 }
01415
01416 if ( !(report.nan_residual_found || report.too_large_residual_found) ) {
01417
01418 for ( int compIdx = 0; compIdx < numComp; ++compIdx )
01419 {
01420 report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells);
01421 }
01422 } else {
01423 report.converged = false;
01424 }
01425
01426 return report;
01427 }
01428
01429
01430
01431
01432
01433 template<typename TypeTag>
01434 void
01435 StandardWell<TypeTag>::
01436 computeWellConnectionDensitesPressures(const WellState& xw,
01437 const std::vector<double>& b_perf,
01438 const std::vector<double>& rsmax_perf,
01439 const std::vector<double>& rvmax_perf,
01440 const std::vector<double>& surf_dens_perf)
01441 {
01442
01443 const int nperf = number_of_perforations_;
01444 const int numComponent = numComponents();
01445 const int np = number_of_phases_;
01446 std::vector<double> perfRates(b_perf.size(),0.0);
01447
01448 for (int perf = 0; perf < nperf; ++perf) {
01449 for (int phase = 0; phase < np; ++phase) {
01450 perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase];
01451 }
01452 if(has_solvent) {
01453 perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf];
01454 }
01455 }
01456
01457 computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
01458
01459 computeConnectionPressureDelta();
01460
01461 }
01462
01463
01464
01465
01466
01467 template<typename TypeTag>
01468 void
01469 StandardWell<TypeTag>::
01470 computeWellConnectionPressures(const Simulator& ebosSimulator,
01471 const WellState& well_state)
01472 {
01473
01474
01475
01476 std::vector<double> b_perf;
01477 std::vector<double> rsmax_perf;
01478 std::vector<double> rvmax_perf;
01479 std::vector<double> surf_dens_perf;
01480 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
01481 computeWellConnectionDensitesPressures(well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
01482 }
01483
01484
01485
01486
01487
01488 template<typename TypeTag>
01489 void
01490 StandardWell<TypeTag>::
01491 solveEqAndUpdateWellState(WellState& well_state)
01492 {
01493
01494
01495 BVectorWell dx_well(1);
01496 invDuneD_.mv(resWell_, dx_well);
01497
01498 updateWellState(dx_well, well_state);
01499 }
01500
01501
01502
01503
01504
01505 template<typename TypeTag>
01506 void
01507 StandardWell<TypeTag>::
01508 calculateExplicitQuantities(const Simulator& ebosSimulator,
01509 const WellState& well_state)
01510 {
01511 computeWellConnectionPressures(ebosSimulator, well_state);
01512 computeAccumWell();
01513 }
01514
01515
01516
01517
01518
01519 template<typename TypeTag>
01520 void
01521 StandardWell<TypeTag>::
01522 computeAccumWell()
01523 {
01524
01525
01526 for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
01527 F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
01528 }
01529 }
01530
01531
01532
01533
01534
01535 template<typename TypeTag>
01536 void
01537 StandardWell<TypeTag>::
01538 apply(const BVector& x, BVector& Ax) const
01539 {
01540 assert( Bx_.size() == duneB_.N() );
01541 assert( invDrw_.size() == invDuneD_.N() );
01542
01543
01544 duneB_.mv(x, Bx_);
01545
01546
01547
01548 BVectorWell& invDBx = invDrw_;
01549 invDuneD_.mv(Bx_, invDBx);
01550
01551
01552 duneC_.mmtv(invDBx,Ax);
01553 }
01554
01555
01556
01557
01558 template<typename TypeTag>
01559 void
01560 StandardWell<TypeTag>::
01561 apply(BVector& r) const
01562 {
01563 assert( invDrw_.size() == invDuneD_.N() );
01564
01565
01566 invDuneD_.mv(resWell_, invDrw_);
01567
01568 duneC_.mmtv(invDrw_, r);
01569 }
01570
01571
01572
01573
01574
01575 template<typename TypeTag>
01576 void
01577 StandardWell<TypeTag>::
01578 recoverSolutionWell(const BVector& x, BVectorWell& xw) const
01579 {
01580 BVectorWell resWell = resWell_;
01581
01582 duneB_.mmv(x, resWell);
01583
01584 invDuneD_.mv(resWell, xw);
01585 }
01586
01587
01588
01589
01590
01591 template<typename TypeTag>
01592 void
01593 StandardWell<TypeTag>::
01594 recoverWellSolutionAndUpdateWellState(const BVector& x,
01595 WellState& well_state) const
01596 {
01597 BVectorWell xw(1);
01598 recoverSolutionWell(x, xw);
01599 updateWellState(xw, well_state);
01600 }
01601
01602
01603
01604
01605
01606 template<typename TypeTag>
01607 void
01608 StandardWell<TypeTag>::
01609 computeWellRatesWithBhp(const Simulator& ebosSimulator,
01610 const EvalWell& bhp,
01611 std::vector<double>& well_flux) const
01612 {
01613 const int np = number_of_phases_;
01614 const int numComp = numComponents();
01615 well_flux.resize(np, 0.0);
01616
01617 const bool allow_cf = crossFlowAllowed(ebosSimulator);
01618
01619 for (int perf = 0; perf < number_of_perforations_; ++perf) {
01620 const int cell_idx = well_cells_[perf];
01621 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
01622
01623 std::vector<EvalWell> cq_s(numComp, 0.0);
01624 std::vector<EvalWell> mob(numComp, 0.0);
01625 getMobility(ebosSimulator, perf, mob);
01626 computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
01627
01628 for(int p = 0; p < np; ++p) {
01629 well_flux[p] += cq_s[p].value();
01630 }
01631 }
01632 }
01633
01634
01635
01636
01637
01638 template<typename TypeTag>
01639 std::vector<double>
01640 StandardWell<TypeTag>::
01641 computeWellPotentialWithTHP(const Simulator& ebosSimulator,
01642 const double initial_bhp,
01643 const std::vector<double>& initial_potential) const
01644 {
01645
01646
01647 const int np = number_of_phases_;
01648
01649 assert( np == int(initial_potential.size()) );
01650
01651 std::vector<double> potentials = initial_potential;
01652 std::vector<double> old_potentials = potentials;
01653
01654 double bhp = initial_bhp;
01655 double old_bhp = bhp;
01656
01657 bool converged = false;
01658 const int max_iteration = 1000;
01659 const double bhp_tolerance = 1000.;
01660
01661 int iteration = 0;
01662
01663 while ( !converged && iteration < max_iteration ) {
01664
01665
01666
01667
01668 bhp = initial_bhp;
01669
01670
01671 const int nwc = well_controls_get_num(well_controls_);
01672
01673 for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
01674 if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
01675 const Opm::PhaseUsage& pu = phaseUsage();
01676
01677 std::vector<double> rates(3, 0.0);
01678 if (active()[ Water ]) {
01679 rates[ Water ] = potentials[pu.phase_pos[ Water ] ];
01680 }
01681 if (active()[ Oil ]) {
01682 rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ];
01683 }
01684 if (active()[ Gas ]) {
01685 rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ];
01686 }
01687
01688 const double bhp_calculated = calculateBhpFromThp(rates, ctrl_index);
01689
01690 if (well_type_ == INJECTOR && bhp_calculated < bhp ) {
01691 bhp = bhp_calculated;
01692 }
01693
01694 if (well_type_ == PRODUCER && bhp_calculated > bhp) {
01695 bhp = bhp_calculated;
01696 }
01697 }
01698 }
01699
01700
01701 if (std::isinf(bhp) || std::isnan(bhp)) {
01702 OPM_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << name());
01703 }
01704
01705 converged = std::abs(old_bhp - bhp) < bhp_tolerance;
01706
01707 computeWellRatesWithBhp(ebosSimulator, bhp, potentials);
01708
01709
01710 for (const double value : potentials) {
01711 if (std::isinf(value) || std::isnan(value)) {
01712 OPM_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << name());
01713 }
01714 }
01715
01716 if (!converged) {
01717 old_bhp = bhp;
01718 for (int p = 0; p < np; ++p) {
01719
01720
01721 const double potential_update_damping_factor = 0.001;
01722 potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p];
01723 old_potentials[p] = potentials[p];
01724 }
01725 }
01726
01727 ++iteration;
01728 }
01729
01730 if (!converged) {
01731 OPM_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << name());
01732 }
01733
01734 return potentials;
01735 }
01736
01737
01738
01739
01740
01741 template<typename TypeTag>
01742 void
01743 StandardWell<TypeTag>::
01744 computeWellPotentials(const Simulator& ebosSimulator,
01745 const WellState& well_state,
01746 std::vector<double>& well_potentials)
01747 {
01748 updatePrimaryVariables(well_state);
01749 computeWellConnectionPressures(ebosSimulator, well_state);
01750
01751
01752
01753 initPrimaryVariablesEvaluation();
01754
01755 const int np = number_of_phases_;
01756 well_potentials.resize(np, 0.0);
01757
01758
01759 const double bhp = mostStrictBhpFromBhpLimits();
01760
01761
01762 if ( !wellHasTHPConstraints() ) {
01763 assert(std::abs(bhp) != std::numeric_limits<double>::max());
01764
01765 computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
01766 } else {
01767
01768
01769 if ( !well_state.isNewWell(index_of_well_) ) {
01770 for (int p = 0; p < np; ++p) {
01771
01772
01773 well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p];
01774 }
01775 } else {
01776
01777 computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
01778 for (double& value : well_potentials) {
01779
01780
01781 const double rate_safety_scaling_factor = 0.00001;
01782 value *= rate_safety_scaling_factor;
01783 }
01784 }
01785
01786 well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials);
01787 }
01788 }
01789
01790
01791
01792
01793
01794 template<typename TypeTag>
01795 void
01796 StandardWell<TypeTag>::
01797 updatePrimaryVariables(const WellState& well_state) const
01798 {
01799 const int np = number_of_phases_;
01800 const int well_index = index_of_well_;
01801 const WellControls* wc = well_controls_;
01802 const double* distr = well_controls_get_current_distr(wc);
01803 const auto pu = phaseUsage();
01804
01805 switch (well_controls_get_current_type(wc)) {
01806 case THP:
01807 case BHP: {
01808 primary_variables_[XvarWell] = 0.0;
01809 if (well_type_ == INJECTOR) {
01810 for (int p = 0; p < np; ++p) {
01811 primary_variables_[XvarWell] += well_state.wellRates()[np*well_index + p] * comp_frac_[p];
01812 }
01813 } else {
01814 for (int p = 0; p < np; ++p) {
01815 primary_variables_[XvarWell] += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
01816 }
01817 }
01818 break;
01819 }
01820 case RESERVOIR_RATE:
01821 case SURFACE_RATE:
01822 primary_variables_[XvarWell] = well_state.bhp()[well_index];
01823 break;
01824 }
01825
01826 double tot_well_rate = 0.0;
01827 for (int p = 0; p < np; ++p) {
01828 tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
01829 }
01830 if(std::abs(tot_well_rate) > 0) {
01831 if (active()[ Water ]) {
01832 primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
01833 }
01834 if (active()[ Gas ]) {
01835 primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
01836 }
01837 if (has_solvent) {
01838 primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / tot_well_rate ;
01839 }
01840 } else {
01841 if (well_type_ == INJECTOR) {
01842
01843 if (active()[Water]) {
01844 if (distr[Water] > 0.0) {
01845 primary_variables_[WFrac] = 1.0;
01846 } else {
01847 primary_variables_[WFrac] = 0.0;
01848 }
01849 }
01850
01851 if (active()[Gas]) {
01852 if (distr[pu.phase_pos[Gas]] > 0.0) {
01853 primary_variables_[GFrac] = 1.0 - wsolvent();
01854 if (has_solvent) {
01855 primary_variables_[SFrac] = wsolvent();
01856 }
01857 } else {
01858 primary_variables_[GFrac] = 0.0;
01859 }
01860 }
01861
01862
01863
01864
01865 } else if (well_type_ == PRODUCER) {
01866
01867 if (active()[Water]) {
01868 primary_variables_[WFrac] = 1.0 / np;
01869 }
01870 if (active()[Gas]) {
01871 primary_variables_[GFrac] = 1.0 / np;
01872 }
01873 } else {
01874 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
01875 }
01876 }
01877 }
01878
01879
01880
01881
01882
01883 template<typename TypeTag>
01884 template<class ValueType>
01885 ValueType
01886 StandardWell<TypeTag>::
01887 calculateBhpFromThp(const std::vector<ValueType>& rates,
01888 const int control_index) const
01889 {
01890
01891
01892
01893
01894
01895
01896 assert(int(rates.size()) == 3);
01897
01898 const ValueType aqua = rates[Water];
01899 const ValueType liquid = rates[Oil];
01900 const ValueType vapour = rates[Gas];
01901
01902 const int vfp = well_controls_iget_vfp(well_controls_, control_index);
01903 const double& thp = well_controls_iget_target(well_controls_, control_index);
01904 const double& alq = well_controls_iget_alq(well_controls_, control_index);
01905
01906
01907 const double rho = perf_densities_[0];
01908
01909 ValueType bhp = 0.;
01910 if (well_type_ == INJECTOR) {
01911 const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
01912
01913 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
01914
01915 bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
01916 }
01917 else if (well_type_ == PRODUCER) {
01918 const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
01919
01920 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
01921
01922 bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
01923 }
01924 else {
01925 OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
01926 }
01927
01928 return bhp;
01929 }
01930
01931
01932
01933
01934
01935 template<typename TypeTag>
01936 double
01937 StandardWell<TypeTag>::
01938 calculateThpFromBhp(const std::vector<double>& rates,
01939 const int control_index,
01940 const double bhp) const
01941 {
01942 assert(int(rates.size()) == 3);
01943
01944 const double aqua = rates[Water];
01945 const double liquid = rates[Oil];
01946 const double vapour = rates[Gas];
01947
01948 const int vfp = well_controls_iget_vfp(well_controls_, control_index);
01949 const double& alq = well_controls_iget_alq(well_controls_, control_index);
01950
01951
01952 const double rho = perf_densities_[0];
01953
01954 double thp = 0.0;
01955 if (well_type_ == INJECTOR) {
01956 const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
01957
01958 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
01959
01960 thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp);
01961 }
01962 else if (well_type_ == PRODUCER) {
01963 const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
01964
01965 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
01966
01967 thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq);
01968 }
01969 else {
01970 OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
01971 }
01972
01973 return thp;
01974 }
01975
01976 template<typename TypeTag>
01977 double
01978 StandardWell<TypeTag>::scalingFactor(const int phaseIdx) const
01979 {
01980 const WellControls* wc = well_controls_;
01981 const double* distr = well_controls_get_current_distr(wc);
01982
01983 if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
01984 if (has_solvent && phaseIdx == contiSolventEqIdx )
01985 OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented");
01986
01987 return distr[phaseIdx];
01988 }
01989 const auto& pu = phaseUsage();
01990 if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
01991 return 1.0;
01992 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
01993 return 1.0;
01994 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
01995 return 0.01;
01996 if (has_solvent && phaseIdx == contiSolventEqIdx )
01997 return 0.01;
01998
01999
02000 assert(false);
02001 return 1.0;
02002 }
02003
02004
02005 }