00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
00022 #define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
00023
00024
00025 namespace Opm
00026 {
00027
00028
00029
00030 namespace wellhelpers {
00031
00032 using ADB = MultisegmentWells::ADB;
00033 using Vector = MultisegmentWells::Vector;
00034
00035 inline
00036 ADB onlyWellDerivs(const ADB& x)
00037 {
00038 Vector val = x.value();
00039 const int nb = x.numBlocks();
00040 if (nb < 2) {
00041 OPM_THROW(std::logic_error, "Called onlyWellDerivs() with argument that has " << nb << " blocks.");
00042 }
00043 std::vector<ADB::M> derivs = { x.derivative()[nb - 2], x.derivative()[nb - 1] };
00044 return ADB::function(std::move(val), std::move(derivs));
00045 }
00046 }
00047
00048
00049
00050
00051
00052 template <class ReservoirResidualQuant, class SolutionState>
00053 void
00054 MultisegmentWells::
00055 extractWellPerfProperties(const SolutionState& ,
00056 const std::vector<ReservoirResidualQuant>& rq,
00057 std::vector<ADB>& mob_perfcells,
00058 std::vector<ADB>& b_perfcells) const
00059 {
00060
00061
00062 if ( !localWellsActive() ) {
00063 mob_perfcells.clear();
00064 b_perfcells.clear();
00065 return;
00066 } else {
00067 const std::vector<int>& well_cells = wellOps().well_cells;
00068 mob_perfcells.resize(num_phases_, ADB::null());
00069 b_perfcells.resize(num_phases_, ADB::null());
00070 for (int phase = 0; phase < num_phases_; ++phase) {
00071 mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
00072 b_perfcells[phase] = subset(rq[phase].b, well_cells);
00073 }
00074 }
00075 }
00076
00077
00078
00079
00080
00081 template <class WellState>
00082 void
00083 MultisegmentWells::
00084 updateWellState(const Vector& dwells,
00085 const double dpmaxrel,
00086 WellState& well_state) const
00087 {
00088 if (!msWells().empty())
00089 {
00090 const int nw = msWells().size();
00091 const int nseg_total = nseg_total_;
00092 const int np = numPhases();
00093
00094
00095 int varstart = 0;
00096 const Vector dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart));
00097 varstart += dsegqs.size();
00098 const Vector dsegp = subset(dwells, Span(nseg_total, 1, varstart));
00099 varstart += dsegp.size();
00100 assert(varstart == dwells.size());
00101
00102
00103
00104
00105
00106 const DataBlock wsr = Eigen::Map<const DataBlock>(dsegqs.data(), np, nseg_total).transpose();
00107 const Vector dwsr = Eigen::Map<const Vector>(wsr.data(), nseg_total * np);
00108 const Vector wsr_old = Eigen::Map<const Vector>(&well_state.segPhaseRates()[0], nseg_total * np);
00109 const Vector sr = wsr_old - dwsr;
00110 std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin());
00111
00112
00113
00114 const Vector segp_old = Eigen::Map<const Vector>(&well_state.segPress()[0], nseg_total, 1);
00115
00116 const Vector dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel);
00117 const Vector segp = segp_old - dsegp_limited;
00118 std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin());
00119
00120
00121
00122
00123
00124 Vector bhp = Vector::Zero(nw);
00125 Vector wr = Vector::Zero(nw * np);
00126
00127
00128 int start_segment = 0;
00129 for (int w = 0; w < nw; ++w) {
00130 bhp[w] = well_state.segPress()[start_segment];
00131
00132 for (int p = 0; p < np; ++p) {
00133 wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
00134 }
00135
00136 const int nseg = msWells()[w]->numberOfSegments();
00137 start_segment += nseg;
00138 }
00139
00140 assert(start_segment == nseg_total);
00141 std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
00142 std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
00143
00144
00145 }
00146 }
00147
00148
00149
00150
00151
00152 template <class SolutionState>
00153 void
00154 MultisegmentWells::
00155 computeWellFlux(const SolutionState& state,
00156 const std::vector<ADB>& mob_perfcells,
00157 const std::vector<ADB>& b_perfcells,
00158 Vector& aliveWells,
00159 std::vector<ADB>& cq_s) const
00160 {
00161 if (msWells().size() == 0) return;
00162
00163 const int np = numPhases();
00164 const int nw = msWells().size();
00165
00166 aliveWells = Vector::Constant(nw, 1.0);
00167
00168 const int nseg = nseg_total_;
00169 const int nperf = nperf_total_;
00170
00171 const Opm::PhaseUsage& pu = fluid_->phaseUsage();
00172
00173 cq_s.resize(np, ADB::null());
00174
00175 {
00176 const Vector& Tw = wellOps().conn_trans_factors;
00177 const std::vector<int>& well_cells = wellOps().well_cells;
00178
00179
00180
00181 const ADB& p_perfcells = subset(state.pressure, well_cells);
00182 const ADB& rs_perfcells = subset(state.rs, well_cells);
00183 const ADB& rv_perfcells = subset(state.rv, well_cells);
00184
00185 const ADB& seg_pressures = state.segp;
00186
00187 const ADB seg_pressures_perf = wellOps().s2p * seg_pressures;
00188
00189
00190 Vector is_multisegment_well(nw);
00191 for (int w = 0; w < nw; ++w) {
00192 is_multisegment_well[w] = double(msWells()[w]->isMultiSegmented());
00193 }
00194
00195 Vector is_multisegment_perf = wellOps().w2p * is_multisegment_well.matrix();
00196 Selector<double> msperf_selector(is_multisegment_perf, Selector<double>::NotEqualZero);
00197
00198
00199 ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
00200 ADB::constant(well_perforation_pressure_diffs_));
00201 const Vector h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, Vector::Zero(nperf));
00202
00203
00204
00205 if ((h_nc.numBlocks() != 0) && (h_nc.numBlocks() != seg_pressures_perf.numBlocks())) {
00206 assert(seg_pressures_perf.numBlocks() == 2);
00207 assert(h_nc.numBlocks() > 2);
00208 h_nc = wellhelpers::onlyWellDerivs(h_nc);
00209 assert(h_nc.numBlocks() == 2);
00210 }
00211
00212 ADB drawdown = (p_perfcells + h_cj - seg_pressures_perf - h_nc);
00213
00214
00215 Vector selectInjectingPerforations = Vector::Zero(nperf);
00216
00217 Vector selectProducingPerforations = Vector::Zero(nperf);
00218 for (int c = 0; c < nperf; ++c){
00219 if (drawdown.value()[c] < 0)
00220 selectInjectingPerforations[c] = 1;
00221 else
00222 selectProducingPerforations[c] = 1;
00223 }
00224
00225
00226
00227 std::vector<ADB> cq_ps(np, ADB::null());
00228 for (int phase = 0; phase < np; ++phase) {
00229 const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
00230 cq_ps[phase] = b_perfcells[phase] * cq_p;
00231 }
00232
00233 if ((*active_)[Oil] && (*active_)[Gas]) {
00234 const int oilpos = pu.phase_pos[Oil];
00235 const int gaspos = pu.phase_pos[Gas];
00236 const ADB cq_psOil = cq_ps[oilpos];
00237 const ADB cq_psGas = cq_ps[gaspos];
00238 cq_ps[gaspos] += rs_perfcells * cq_psOil;
00239 cq_ps[oilpos] += rv_perfcells * cq_psGas;
00240 }
00241
00242
00243 ADB total_mob = mob_perfcells[0];
00244 for (int phase = 1; phase < np; ++phase) {
00245 total_mob += mob_perfcells[phase];
00246 }
00247
00248
00249 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 std::vector<ADB> wbq(np, ADB::null());
00268 ADB wbqt = ADB::constant(Vector::Zero(nseg));
00269
00270 const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
00271
00272 for (int phase = 0; phase < np; ++phase) {
00273 const ADB& q_ps = wellOps().p2s * cq_ps[phase];
00274 const ADB& q_s = subset(state.segqs, Span(nseg, 1, phase * nseg));
00275 Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
00276
00277 const int pos = pu.phase_pos[phase];
00278
00279
00280 wbq[phase] = (wellOps().w2s * ADB::constant(compi.col(pos)) * injectingPhase_selector.select(q_s, ADB::constant(Vector::Zero(nseg)))) - q_ps;
00281
00282
00283
00284
00285 wbqt += wbq[phase];
00286 }
00287
00288
00289
00290
00291 {
00292 int topseg = 0;
00293 for (int w = 0; w < nw; ++w) {
00294 if (wbqt.value()[topseg] == 0.0) {
00295 aliveWells[w] = 0.0;
00296 }
00297 topseg += msWells()[w]->numberOfSegments();
00298 }
00299 }
00300
00301
00302
00303
00304
00305 std::vector<ADB> cmix_s(np, ADB::null());
00306 Selector<double> aliveWells_selector(aliveWells, Selector<double>::NotEqualZero);
00307 for (int phase = 0; phase < np; ++phase) {
00308 const int pos = pu.phase_pos[phase];
00309 const ADB phase_fraction = wellOps().topseg2w * (wbq[phase] / wbqt);
00310 cmix_s[phase] = wellOps().w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos)));
00311 }
00312
00313
00314 ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
00315 const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
00316
00317 for (int phase = 0; phase < np; ++phase) {
00318 ADB tmp = cmix_s[phase];
00319 if (phase == Oil && (*active_)[Gas]) {
00320 const int gaspos = pu.phase_pos[Gas];
00321 tmp = (tmp - rv_perfcells * cmix_s[gaspos]) / d;
00322 }
00323 if (phase == Gas && (*active_)[Oil]) {
00324 const int oilpos = pu.phase_pos[Oil];
00325 tmp = (tmp - rs_perfcells * cmix_s[oilpos]) / d;
00326 }
00327 volumeRatio += tmp / b_perfcells[phase];
00328 }
00329
00330
00331 ADB cqt_is = cqt_i/volumeRatio;
00332
00333
00334 for (int phase = 0; phase < np; ++phase) {
00335 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
00336 }
00337 }
00338 }
00339
00340
00341
00342
00343
00344 template <class SolutionState, class WellState>
00345 void
00346 MultisegmentWells::
00347 updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
00348 const SolutionState& state,
00349 WellState& xw) const
00350 {
00351 if ( !localWellsActive() ) {
00352 return;
00353 }
00354
00355 const int np = numPhases();
00356 const int nw = numWells();
00357
00358 Vector cq = superset(cq_s[0].value(), Span(nperf_total_, np, 0), nperf_total_ * np);
00359 for (int phase = 1; phase < np; ++phase) {
00360 cq += superset(cq_s[phase].value(), Span(nperf_total_, np, phase), nperf_total_ * np);
00361 }
00362 xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total_ * np);
00363
00364
00365
00366
00367 xw.perfPress().resize(nperf_total_, -1.e100);
00368
00369 const Vector& cdp = well_perforation_pressure_diffs_;
00370 int start_segment = 0;
00371 int start_perforation = 0;
00372 for (int i = 0; i < nw; ++i) {
00373 WellMultiSegmentConstPtr well = wells_multisegment_[i];
00374 const int nperf = well->numberOfPerforations();
00375 const int nseg = well->numberOfSegments();
00376 if (well->isMultiSegmented()) {
00377 start_segment += nseg;
00378 start_perforation += nperf;
00379 continue;
00380 }
00381 const Vector cdp_well = subset(cdp, Span(nperf, 1, start_perforation));
00382 const ADB segp = subset(state.segp, Span(nseg, 1, start_segment));
00383 const Vector perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well;
00384 std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]);
00385
00386 start_segment += nseg;
00387 start_perforation += nperf;
00388 }
00389 assert(start_segment == nseg_total_);
00390 assert(start_perforation == nperf_total_);
00391 }
00392
00393
00394
00395
00396
00397 template <class SolutionState>
00398 void
00399 MultisegmentWells::
00400 computeSegmentFluidProperties(const SolutionState& state)
00401 {
00402 const int np = numPhases();
00403 const int nw = msWells().size();
00404 const int nseg_total = nseg_total_;
00405
00406 if ( !wellOps().has_multisegment_wells ){
00407
00408
00409 well_segment_densities_ = ADB::constant(Vector::Zero(nseg_total));
00410 segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
00411 segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
00412 for (int phase = 0; phase < np; ++phase) {
00413 segment_comp_surf_volume_current_[phase] = ADB::constant(Vector::Zero(nseg_total));
00414 segmentCompSurfVolumeInitial()[phase] = Vector::Zero(nseg_total);
00415 }
00416 return;
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428 std::vector<int> segment_cells;
00429 segment_cells.reserve(nseg_total);
00430 for (int w = 0; w < nw; ++w) {
00431 const std::vector<int>& segment_cells_well = msWells()[w]->segmentCells();
00432 segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end());
00433 }
00434 assert(int(segment_cells.size()) == nseg_total);
00435
00436 const ADB segment_temp = subset(state.temperature, segment_cells);
00437
00438
00439 const ADB& segment_press = state.segp;
00440
00441
00442 std::vector<PhasePresence> segment_cond(nseg_total);
00443 for (int s = 0; s < nseg_total; ++s) {
00444 segment_cond[s] = (*phase_condition_)[segment_cells[s]];
00445 }
00446 std::vector<ADB> b_seg(np, ADB::null());
00447
00448 std::vector<ADB> mu_seg(np, ADB::null());
00449 ADB rsmax_seg = ADB::null();
00450 ADB rvmax_seg = ADB::null();
00451 const PhaseUsage& pu = fluid_->phaseUsage();
00452 if (pu.phase_used[Water]) {
00453 b_seg[pu.phase_pos[Water]] = fluid_->bWat(segment_press, segment_temp, segment_cells);
00454 mu_seg[pu.phase_pos[Water]] = fluid_->muWat(segment_press, segment_temp, segment_cells);
00455 }
00456 assert((*active_)[Oil]);
00457 const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
00458 if (pu.phase_used[Oil]) {
00459 const ADB segment_rs = subset(state.rs, segment_cells);
00460 b_seg[pu.phase_pos[Oil]] = fluid_->bOil(segment_press, segment_temp, segment_rs,
00461 segment_cond, segment_cells);
00462
00463 rsmax_seg = fluid_->rsSat(segment_press, segment_so, segment_cells);
00464 mu_seg[pu.phase_pos[Oil]] = fluid_->muOil(segment_press, segment_temp, segment_rs,
00465 segment_cond, segment_cells);
00466 }
00467 assert((*active_)[Gas]);
00468 if (pu.phase_used[Gas]) {
00469 const ADB segment_rv = subset(state.rv, segment_cells);
00470 b_seg[pu.phase_pos[Gas]] = fluid_->bGas(segment_press, segment_temp, segment_rv,
00471 segment_cond, segment_cells);
00472
00473 rvmax_seg = fluid_->rvSat(segment_press, segment_so, segment_cells);
00474 mu_seg[pu.phase_pos[Gas]] = fluid_->muGas(segment_press, segment_temp, segment_rv,
00475 segment_cond, segment_cells);
00476 }
00477
00478
00479 ADB tot_surface_rate = ADB::constant(Vector::Zero(nseg_total));
00480 std::vector<ADB> segqs(np, ADB::null());
00481 for (int phase = 0; phase < np; ++phase) {
00482 segqs[phase] = subset(state.segqs, Span(nseg_total, 1, phase * nseg_total));
00483 tot_surface_rate += segqs[phase];
00484 }
00485
00486
00487 std::vector<std::vector<double>> comp_frac(np, std::vector<double>(nseg_total, 0.0));
00488 int start_segment = 0;
00489 for (int w = 0; w < nw; ++w) {
00490 WellMultiSegmentConstPtr well = msWells()[w];
00491 const int nseg = well->numberOfSegments();
00492 const std::vector<double>& comp_frac_well = well->compFrac();
00493 for (int phase = 0; phase < np; ++phase) {
00494 for (int s = 0; s < nseg; ++s) {
00495 comp_frac[phase][s + start_segment] = comp_frac_well[phase];
00496 }
00497 }
00498 start_segment += nseg;
00499 }
00500 assert(start_segment == nseg_total);
00501
00502
00503
00504 std::vector<ADB> mix(np, ADB::null());
00505 for (int phase = 0; phase < np; ++phase) {
00506
00507
00508 mix[phase] = ADB::constant(Eigen::Map<Vector>(comp_frac[phase].data(), nseg_total));
00509 }
00510
00511 Selector<double> non_zero_tot_rate(tot_surface_rate.value(), Selector<double>::NotEqualZero);
00512 for (int phase = 0; phase < np; ++phase) {
00513 mix[phase] = non_zero_tot_rate.select(segqs[phase] / tot_surface_rate, mix[phase]);
00514 }
00515
00516
00517 ADB rs = ADB::constant(Vector::Zero(nseg_total));
00518 ADB rv = rs;
00519 const int gaspos = pu.phase_pos[Gas];
00520 const int oilpos = pu.phase_pos[Oil];
00521 Selector<double> non_zero_mix_oilpos(mix[oilpos].value(), Selector<double>::GreaterZero);
00522 Selector<double> non_zero_mix_gaspos(mix[gaspos].value(), Selector<double>::GreaterZero);
00523
00524
00525 ADB big_values = ADB::constant(Vector::Constant(nseg_total, 1.e100));
00526 ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
00527 ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
00528 if ((*active_)[Oil]) {
00529 Vector selectorUnderRsmax = Vector::Zero(nseg_total);
00530 Vector selectorAboveRsmax = Vector::Zero(nseg_total);
00531 for (int s = 0; s < nseg_total; ++s) {
00532 if (mix_gas_oil.value()[s] > rsmax_seg.value()[s]) {
00533 selectorAboveRsmax[s] = 1.0;
00534 } else {
00535 selectorUnderRsmax[s] = 1.0;
00536 }
00537 }
00538 rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
00539 }
00540 if ((*active_)[Gas]) {
00541 Vector selectorUnderRvmax = Vector::Zero(nseg_total);
00542 Vector selectorAboveRvmax = Vector::Zero(nseg_total);
00543 for (int s = 0; s < nseg_total; ++s) {
00544 if (mix_oil_gas.value()[s] > rvmax_seg.value()[s]) {
00545 selectorAboveRvmax[s] = 1.0;
00546 } else {
00547 selectorUnderRvmax[s] = 1.0;
00548 }
00549 }
00550 rv = non_zero_mix_gaspos.select(selectorAboveRvmax * rvmax_seg + selectorUnderRvmax * mix_oil_gas, rv);
00551 }
00552
00553
00554 std::vector<ADB> x(np, ADB::null());
00555 for (int phase = 0; phase < np; ++phase) {
00556 x[phase] = mix[phase];
00557 }
00558 if ((*active_)[Gas] && (*active_)[Oil]) {
00559 x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
00560 x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
00561 }
00562
00563
00564 ADB volrat = ADB::constant(Vector::Zero(nseg_total));
00565 for (int phase = 0; phase < np; ++phase) {
00566 volrat += x[phase] / b_seg[phase];
00567 }
00568
00569
00570 ADB dens = ADB::constant(Vector::Zero(nseg_total));
00571 for (int phase = 0; phase < np; ++phase) {
00572 const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
00573 dens += surface_density * mix[phase];
00574 }
00575 well_segment_densities_ = dens / volrat;
00576
00577
00578 assert(np == int(segment_comp_surf_volume_current_.size()));
00579 const ADB segment_surface_volume = segvdt_ / volrat;
00580 for (int phase = 0; phase < np; ++phase) {
00581 segment_comp_surf_volume_current_[phase] = segment_surface_volume * mix[phase];
00582 }
00583
00584
00585 segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
00586 for (int phase = 0; phase < np; ++phase) {
00587
00588 const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
00589 segment_mass_flow_rates_ += surface_density * segqs[phase];
00590 }
00591
00592
00593 segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
00594 for (int phase = 0; phase < np; ++phase) {
00595 segment_viscosities_ += x[phase] * mu_seg[phase];
00596 }
00597 }
00598
00599
00600
00601
00602
00603 template <class SolutionState>
00604 void
00605 MultisegmentWells::
00606 addWellFluxEq(const std::vector<ADB>& cq_s,
00607 const SolutionState& state,
00608 LinearisedBlackoilResidual& residual)
00609 {
00610 if ( !localWellsActive() ) {
00611 return;
00612 }
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 const int np = numPhases();
00625 const int nseg_total = nseg_total_;
00626
00627 ADB segqs = state.segqs;
00628
00629 std::vector<ADB> segment_volume_change_dt(np, ADB::null());
00630 for (int phase = 0; phase < np; ++phase) {
00631 if ( wellOps().has_multisegment_wells ) {
00632
00633 segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
00634 segmentCompSurfVolumeInitial()[phase];
00635
00636
00637
00638 if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
00639 assert(segment_volume_change_dt[phase].numBlocks() > 2);
00640 assert(segqs.numBlocks() == 2);
00641 segment_volume_change_dt[phase] = wellhelpers::onlyWellDerivs(segment_volume_change_dt[phase]);
00642 assert(segment_volume_change_dt[phase].numBlocks() == 2);
00643 }
00644
00645 const ADB cq_s_seg = wellOps().p2s * cq_s[phase];
00646 const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total));
00647 segqs -= superset(cq_s_seg + wellOps().s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
00648 Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
00649 } else {
00650 segqs -= superset(wellOps().p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
00651 }
00652 }
00653
00654 residual.well_flux_eq = segqs;
00655 }
00656
00657
00658
00659
00660
00661 template <class SolutionState, class WellState>
00662 void
00663 MultisegmentWells::
00664 addWellControlEq(const SolutionState& state,
00665 const WellState& xw,
00666 const Vector& aliveWells,
00667 LinearisedBlackoilResidual& residual)
00668 {
00669
00670
00671
00672 if( msWells().empty() ) return;
00673
00674 const int np = numPhases();
00675 const int nw = msWells().size();
00676 const int nseg_total = nseg_total_;
00677
00678 ADB aqua = ADB::constant(Vector::Zero(nseg_total));
00679 ADB liquid = ADB::constant(Vector::Zero(nseg_total));
00680 ADB vapour = ADB::constant(Vector::Zero(nseg_total));
00681
00682 if ((*active_)[Water]) {
00683 aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
00684 }
00685 if ((*active_)[Oil]) {
00686 liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
00687 }
00688 if ((*active_)[Gas]) {
00689 vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
00690 }
00691
00692
00693
00694
00695 Vector rho_v = Vector::Zero(nw);
00696 Vector vfp_ref_depth_v = Vector::Zero(nw);
00697
00698
00699 Vector bhp_targets = Vector::Zero(nw);
00700 Vector rate_targets = Vector::Zero(nw);
00701 Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
00702
00703
00704
00705 std::vector<int> bhp_well_elems;
00706 std::vector<int> rate_well_elems;
00707
00708 std::vector<int> bhp_top_elems;
00709 std::vector<int> rate_top_elems;
00710 std::vector<int> rate_top_phase_elems;
00711 std::vector<int> others_elems;
00712
00713
00714
00715 int start_segment = 0;
00716 for (int w = 0; w < nw; ++w) {
00717 const struct WellControls* wc = msWells()[w]->wellControls();
00718
00719
00720
00721
00722 const int current = xw.currentControls()[w];
00723
00724 const int nseg = msWells()[w]->numberOfSegments();
00725
00726 switch (well_controls_iget_type(wc, current)) {
00727 case BHP:
00728 {
00729 bhp_well_elems.push_back(w);
00730 bhp_top_elems.push_back(start_segment);
00731 bhp_targets(w) = well_controls_iget_target(wc, current);
00732 rate_targets(w) = -1e100;
00733 for (int p = 0; p < np; ++p) {
00734 rate_top_phase_elems.push_back(np * start_segment + p);
00735 }
00736 }
00737 break;
00738
00739 case THP:
00740 {
00741 OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
00742 }
00743 break;
00744
00745 case RESERVOIR_RATE:
00746 case SURFACE_RATE:
00747 {
00748 rate_well_elems.push_back(w);
00749 rate_top_elems.push_back(start_segment);
00750 for (int p = 0; p < np; ++p) {
00751 rate_top_phase_elems.push_back(np * start_segment + p);
00752 }
00753
00754
00755
00756
00757 const double* const distr =
00758 well_controls_iget_distr(wc, current);
00759
00760 for (int p = 0; p < np; ++p) {
00761 rate_distr.insert(w, p*nw + w) = distr[p];
00762 }
00763
00764 bhp_targets(w) = -1.0e100;
00765 rate_targets(w) = well_controls_iget_target(wc, current);
00766 }
00767 break;
00768
00769 }
00770
00771 for (int i = 1; i < nseg; ++i) {
00772 others_elems.push_back(i + start_segment);
00773 }
00774 start_segment += nseg;
00775 }
00776
00777
00778
00779 const ADB bhp_residual = subset(state.segp, bhp_top_elems) - subset(bhp_targets, bhp_well_elems);
00780 const ADB rate_residual = subset(rate_distr * subset(state.segqs, rate_top_phase_elems) - rate_targets, rate_well_elems);
00781
00782 ADB others_residual = ADB::constant(Vector::Zero(nseg_total));
00783
00784 if ( wellOps().has_multisegment_wells ) {
00785
00786
00787 ADB wspd = (state.segp.numBlocks() == 2)
00788 ? wellhelpers::onlyWellDerivs(well_segment_pressures_delta_)
00789 : well_segment_pressures_delta_;
00790
00791 others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp + wspd);
00792 } else {
00793 others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp);
00794 }
00795
00796
00797
00798 ADB well_eq_topsegment = subset(superset(bhp_residual, bhp_top_elems, nseg_total) +
00799 superset(rate_residual, rate_top_elems, nseg_total), top_well_segments_);
00800
00801
00802
00803
00804 Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
00805 for (int w = 0; w < nw; ++w) {
00806 for (int phase = 0; phase < np; ++phase) {
00807 rate_summer.insert(w, phase*nw + w) = 1.0;
00808 }
00809 }
00810 Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
00811
00812
00813
00814
00815 well_eq_topsegment = alive_selector.select(well_eq_topsegment, rate_summer * subset(state.segqs, rate_top_phase_elems));
00816
00817
00818
00819
00820 residual.well_eq = superset(well_eq_topsegment, top_well_segments_, nseg_total) +
00821 others_residual;
00822 }
00823
00824
00825
00826
00827
00828 template <class WellState>
00829 void
00830 MultisegmentWells::
00831 updateWellControls(WellState& xw) const
00832 {
00833 wellhelpers::WellSwitchingLogger logger;
00834
00835 if( msWells().empty() ) return ;
00836
00837
00838
00839 const int np = numPhases();
00840 const int nw = msWells().size();
00841 for (int w = 0; w < nw; ++w) {
00842 const WellControls* wc = msWells()[w]->wellControls();
00843
00844
00845
00846 int current = xw.currentControls()[w];
00847
00848
00849
00850 const int nwc = well_controls_get_num(wc);
00851 int ctrl_index = 0;
00852 for (; ctrl_index < nwc; ++ctrl_index) {
00853 if (ctrl_index == current) {
00854
00855
00856
00857 continue;
00858 }
00859 if (wellhelpers::constraintBroken(
00860 xw.bhp(), xw.thp(), xw.wellRates(),
00861 w, np, msWells()[w]->wellType(), wc, ctrl_index)) {
00862
00863 break;
00864 }
00865 }
00866
00867 if (ctrl_index != nwc) {
00868
00869
00870 logger.wellSwitched(msWells()[w]->name(),
00871 well_controls_iget_type(wc, current),
00872 well_controls_iget_type(wc, ctrl_index));
00873 xw.currentControls()[w] = ctrl_index;
00874 current = xw.currentControls()[w];
00875 }
00876
00877
00878
00879
00880
00881
00882 const double target = well_controls_iget_target(wc, current);
00883 const double* distr = well_controls_iget_distr(wc, current);
00884 switch (well_controls_iget_type(wc, current)) {
00885 case BHP:
00886 xw.bhp()[w] = target;
00887 xw.segPress()[top_well_segments_[w]] = target;
00888 break;
00889
00890 case THP: {
00891 OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
00892 }
00893
00894 case RESERVOIR_RATE:
00895
00896
00897
00898
00899 break;
00900
00901 case SURFACE_RATE:
00902 for (int phase = 0; phase < np; ++phase) {
00903 if (distr[phase] > 0.0) {
00904 xw.wellRates()[np * w + phase] = target * distr[phase];
00905
00906
00907 xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
00908 }
00909 }
00910 break;
00911 }
00912
00913 if (wellCollection()->groupControlActive()) {
00914
00915 WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
00916
00917
00918 if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
00919
00920 well_node.setIndividualControl(false);
00921 } else {
00922
00923 well_node.setIndividualControl(true);
00924 }
00925 }
00926
00927 }
00928
00929 }
00930
00931
00932
00933
00934
00935
00936
00937 template <class SolutionState, class WellState>
00938 void
00939 MultisegmentWells::
00940 computeWellConnectionPressures(const SolutionState& state,
00941 const WellState& xw,
00942 const std::vector<ADB>& kr_adb,
00943 const std::vector<ADB>& fluid_density)
00944 {
00945 if( ! wellsActive() ) return ;
00946
00947 using namespace Opm::AutoDiffGrid;
00948
00949
00950
00951 const int nperf_total = nperf_total_;
00952 const int nw = numWells();
00953
00954 const std::vector<int>& well_cells = wellOps().well_cells;
00955
00956 well_perforation_densities_ = Vector::Zero(nperf_total);
00957
00958 const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf_total);
00959
00960 Vector avg_press = perf_press * 0.0;
00961
00962
00963
00964
00965 int start_segment = 0;
00966 for (int w = 0; w < nw; ++w) {
00967 const int nseg = wells_multisegment_[w]->numberOfSegments();
00968 if (wells_multisegment_[w]->isMultiSegmented()) {
00969
00970 start_segment += nseg;
00971 continue;
00972 }
00973
00974 std::string well_name(wells_multisegment_[w]->name());
00975 typedef typename WellState::SegmentedWellMapType::const_iterator const_iterator;
00976 const_iterator it_well = xw.segmentedWellMap().find(well_name);
00977 assert(it_well != xw.segmentedWellMap().end());
00978
00979 const int start_perforation = (*it_well).second.start_perforation;
00980 const int end_perforation = start_perforation + (*it_well).second.number_of_perforations;
00981 for (int perf = start_perforation; perf < end_perforation; ++perf) {
00982 const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1];
00983 const double p_avg = (perf_press[perf] + p_above)/2;
00984 avg_press[perf] = p_avg;
00985 }
00986 start_segment += nseg;
00987 }
00988 assert(start_segment == xw.numSegments());
00989
00990
00991 const ADB perf_temp = subset(state.temperature, well_cells);
00992
00993
00994
00995
00996 const ADB avg_press_ad = ADB::constant(avg_press);
00997 std::vector<PhasePresence> perf_cond(nperf_total);
00998 for (int perf = 0; perf < nperf_total; ++perf) {
00999 perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
01000 }
01001 const PhaseUsage& pu = fluid_->phaseUsage();
01002 DataBlock b(nperf_total, pu.num_phases);
01003 std::vector<double> rsmax_perf(nperf_total, 0.0);
01004 std::vector<double> rvmax_perf(nperf_total, 0.0);
01005 if (pu.phase_used[BlackoilPhases::Aqua]) {
01006 const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
01007 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
01008 }
01009 assert((*active_)[Oil]);
01010 const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
01011 if (pu.phase_used[BlackoilPhases::Liquid]) {
01012 const ADB perf_rs = subset(state.rs, well_cells);
01013 const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
01014 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
01015 const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
01016 rsmax_perf.assign(rssat.data(), rssat.data() + nperf_total);
01017 }
01018 if (pu.phase_used[BlackoilPhases::Vapour]) {
01019 const ADB perf_rv = subset(state.rv, well_cells);
01020 const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
01021 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
01022 const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
01023 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf_total);
01024 }
01025
01026 std::vector<double> b_perf(b.data(), b.data() + nperf_total * pu.num_phases);
01027
01028 const Vector& perfcelldepth = perf_cell_depth_;
01029 std::vector<double> perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total);
01030
01031
01032
01033
01034 Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases);
01035 for (int phase = 1; phase < pu.num_phases; ++phase) {
01036 rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases);
01037 }
01038 std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases);
01039
01040
01041 std::vector<double> cd =
01042 WellDensitySegmented::computeConnectionDensities(
01043 wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
01044 b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
01045
01046
01047 std::vector<double> cdp =
01048 WellDensitySegmented::computeConnectionPressureDelta(
01049 wells(), perf_cell_depth, cd, gravity_);
01050
01051
01052 well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf_total);
01053 well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf_total);
01054
01055 if ( !wellOps().has_multisegment_wells ) {
01056 well_perforation_cell_densities_ = Vector::Zero(nperf_total);
01057 well_perforation_cell_pressure_diffs_ = Vector::Zero(nperf_total);
01058 return;
01059 }
01060
01061
01062
01063
01064 size_t temp_size = kr_adb.size();
01065 std::vector<Vector> perf_kr;
01066 for(size_t i = 0; i < temp_size; ++i) {
01067
01068 const Vector kr_phase = (subset(kr_adb[i], well_cells)).value();
01069 perf_kr.push_back(kr_phase);
01070 }
01071
01072
01073
01074
01075
01076 for (int i = 0; i < nperf_total; ++i) {
01077 double sum_kr = 0.;
01078 int np = perf_kr.size();
01079 for (int p = 0; p < np; ++p) {
01080 sum_kr += perf_kr[p][i];
01081 }
01082
01083 for (int p = 0; p < np; ++p) {
01084 perf_kr[p][i] /= sum_kr;
01085 }
01086 }
01087
01088 Vector rho_avg_perf = Vector::Constant(nperf_total, 0.0);
01089
01090 for (int phaseIdx = 0; phaseIdx < fluid_->numPhases(); ++phaseIdx) {
01091
01092
01093 const Vector rho_perf = subset(fluid_density[phaseIdx], well_cells).value();
01094
01095 rho_avg_perf += rho_perf * perf_kr[phaseIdx];
01096 }
01097
01098 well_perforation_cell_densities_ = Eigen::Map<const Vector>(rho_avg_perf.data(), nperf_total);
01099
01100 well_perforation_cell_pressure_diffs_ = gravity_ * well_perforation_cell_densities_ * perf_cell_depth_diffs_;
01101 }
01102
01103
01104
01105
01106
01107 template <class SolutionState>
01108 void
01109 MultisegmentWells::
01110 variableStateExtractWellsVars(const std::vector<int>& indices,
01111 std::vector<ADB>& vars,
01112 SolutionState& state) const
01113 {
01114
01115
01116
01117
01118 state.segqs = std::move(vars[indices[Qs]]);
01119
01120
01121 state.segp = std::move(vars[indices[Bhp]]);
01122
01123
01124
01125
01126 const int np = num_phases_;
01127 const int ns = nseg_total_;
01128 const int nw = numWells();
01129 state.qs = ADB::constant(Vector::Zero(np * nw));
01130 for (int phase = 0; phase < np; ++phase) {
01131
01132 ADB segqs_phase = subset(state.segqs, Span(ns, 1, ns*phase));
01133
01134 ADB wellqs_phase = subset(segqs_phase, topWellSegments());
01135
01136 state.qs += superset(wellqs_phase, Span(nw, 1, nw * phase), nw * np);
01137 }
01138 state.bhp = subset(state.segp, topWellSegments());
01139 }
01140
01141
01142
01143
01144
01145 template <class WellState>
01146 void
01147 MultisegmentWells::
01148 variableWellStateInitials(const WellState& xw,
01149 std::vector<Vector>& vars0) const
01150 {
01151
01152 if ( wells_multisegment_.size() > 0 )
01153 {
01154
01155 const int nseg = xw.numSegments();
01156 const int np = xw.numPhases();
01157
01158
01159 const DataBlock segrates = Eigen::Map<const DataBlock>(& xw.segPhaseRates()[0], nseg, np).transpose();
01160
01161 const Vector segqs = Eigen::Map<const Vector>(segrates.data(), nseg * np);
01162 vars0.push_back(segqs);
01163
01164
01165 const Vector segp = Eigen::Map<const Vector>(& xw.segPress()[0], xw.segPress().size());
01166 vars0.push_back(segp);
01167 }
01168 else
01169 {
01170
01171 vars0.push_back(Vector());
01172 vars0.push_back(Vector());
01173 }
01174 }
01175
01176
01177 }
01178 #endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED