21 #ifndef OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED 22 #define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED 30 namespace wellhelpers {
32 using ADB = MultisegmentWells::ADB;
33 using Vector = MultisegmentWells::Vector;
36 ADB onlyWellDerivs(
const ADB& x)
38 Vector val = x.value();
39 const int nb = x.numBlocks();
41 OPM_THROW(std::logic_error,
"Called onlyWellDerivs() with argument that has " << nb <<
" blocks.");
43 std::vector<ADB::M> derivs = { x.derivative()[nb - 2], x.derivative()[nb - 1] };
52 template <
class ReservoirRes
idualQuant,
class SolutionState>
55 extractWellPerfProperties(
const SolutionState& ,
56 const std::vector<ReservoirResidualQuant>& rq,
57 std::vector<ADB>& mob_perfcells,
58 std::vector<ADB>& b_perfcells)
const 62 if ( !localWellsActive() ) {
63 mob_perfcells.clear();
67 const std::vector<int>& well_cells = wellOps().well_cells;
68 mob_perfcells.resize(num_phases_,
ADB::null());
69 b_perfcells.resize(num_phases_,
ADB::null());
70 for (
int phase = 0; phase < num_phases_; ++phase) {
71 mob_perfcells[phase] =
subset(rq[phase].mob, well_cells);
72 b_perfcells[phase] =
subset(rq[phase].b, well_cells);
81 template <
class WellState>
84 updateWellState(
const Vector& dwells,
85 const double dpmaxrel,
86 WellState& well_state)
const 88 if (!msWells().empty())
90 const int nw = msWells().size();
91 const int nseg_total = nseg_total_;
92 const int np = numPhases();
96 const Vector dsegqs =
subset(dwells, Span(np * nseg_total, 1, varstart));
97 varstart += dsegqs.size();
98 const Vector dsegp =
subset(dwells, Span(nseg_total, 1, varstart));
99 varstart += dsegp.size();
100 assert(varstart == dwells.size());
106 const DataBlock wsr = Eigen::Map<const DataBlock>(dsegqs.data(), np, nseg_total).transpose();
107 const Vector dwsr = Eigen::Map<const Vector>(wsr.data(), nseg_total * np);
108 const Vector wsr_old = Eigen::Map<const Vector>(&well_state.segPhaseRates()[0], nseg_total * np);
109 const Vector sr = wsr_old - dwsr;
110 std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin());
114 const Vector segp_old = Eigen::Map<const Vector>(&well_state.segPress()[0], nseg_total, 1);
116 const Vector dsegp_limited =
sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel);
117 const Vector segp = segp_old - dsegp_limited;
118 std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin());
124 Vector bhp = Vector::Zero(nw);
125 Vector wr = Vector::Zero(nw * np);
128 int start_segment = 0;
129 for (
int w = 0; w < nw; ++w) {
130 bhp[w] = well_state.segPress()[start_segment];
132 for (
int p = 0; p < np; ++p) {
133 wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
136 const int nseg = msWells()[w]->numberOfSegments();
137 start_segment += nseg;
140 assert(start_segment == nseg_total);
141 std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
142 std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
152 template <
class SolutionState>
155 computeWellFlux(
const SolutionState& state,
156 const std::vector<ADB>& mob_perfcells,
157 const std::vector<ADB>& b_perfcells,
159 std::vector<ADB>& cq_s)
const 161 if (msWells().size() == 0)
return;
163 const int np = numPhases();
164 const int nw = msWells().size();
166 aliveWells = Vector::Constant(nw, 1.0);
168 const int nseg = nseg_total_;
169 const int nperf = nperf_total_;
171 const Opm::PhaseUsage& pu = fluid_->
phaseUsage();
176 const Vector& Tw = wellOps().conn_trans_factors;
177 const std::vector<int>& well_cells = wellOps().well_cells;
181 const ADB& p_perfcells =
subset(state.pressure, well_cells);
182 const ADB& rs_perfcells =
subset(state.rs, well_cells);
183 const ADB& rv_perfcells =
subset(state.rv, well_cells);
185 const ADB& seg_pressures = state.segp;
187 const ADB seg_pressures_perf = wellOps().s2p * seg_pressures;
190 Vector is_multisegment_well(nw);
191 for (
int w = 0; w < nw; ++w) {
192 is_multisegment_well[w] = double(msWells()[w]->isMultiSegmented());
195 Vector is_multisegment_perf = wellOps().w2p * is_multisegment_well.matrix();
196 Selector<double> msperf_selector(is_multisegment_perf, Selector<double>::NotEqualZero);
199 ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
201 const Vector h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, Vector::Zero(nperf));
205 if ((h_nc.numBlocks() != 0) && (h_nc.numBlocks() != seg_pressures_perf.numBlocks())) {
206 assert(seg_pressures_perf.numBlocks() == 2);
207 assert(h_nc.numBlocks() > 2);
208 h_nc = wellhelpers::onlyWellDerivs(h_nc);
209 assert(h_nc.numBlocks() == 2);
212 ADB drawdown = (p_perfcells + h_cj - seg_pressures_perf - h_nc);
215 Vector selectInjectingPerforations = Vector::Zero(nperf);
217 Vector selectProducingPerforations = Vector::Zero(nperf);
218 for (
int c = 0; c < nperf; ++c){
219 if (drawdown.value()[c] < 0)
220 selectInjectingPerforations[c] = 1;
222 selectProducingPerforations[c] = 1;
228 for (
int phase = 0; phase < np; ++phase) {
229 const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
230 cq_ps[phase] = b_perfcells[phase] * cq_p;
233 if ((*active_)[Oil] && (*active_)[Gas]) {
234 const int oilpos = pu.phase_pos[Oil];
235 const int gaspos = pu.phase_pos[Gas];
236 const ADB cq_psOil = cq_ps[oilpos];
237 const ADB cq_psGas = cq_ps[gaspos];
238 cq_ps[gaspos] += rs_perfcells * cq_psOil;
239 cq_ps[oilpos] += rv_perfcells * cq_psGas;
243 ADB total_mob = mob_perfcells[0];
244 for (
int phase = 1; phase < np; ++phase) {
245 total_mob += mob_perfcells[phase];
249 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
270 const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
272 for (
int phase = 0; phase < np; ++phase) {
273 const ADB& q_ps = wellOps().p2s * cq_ps[phase];
274 const ADB& q_s =
subset(state.segqs, Span(nseg, 1, phase * nseg));
275 Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
277 const int pos = pu.phase_pos[phase];
280 wbq[phase] = (wellOps().w2s *
ADB::constant(compi.col(pos)) * injectingPhase_selector.select(q_s,
ADB::constant(Vector::Zero(nseg)))) - q_ps;
293 for (
int w = 0; w < nw; ++w) {
294 if (wbqt.value()[topseg] == 0.0) {
297 topseg += msWells()[w]->numberOfSegments();
305 std::vector<ADB> cmix_s(np,
ADB::null());
306 Selector<double> aliveWells_selector(aliveWells, Selector<double>::NotEqualZero);
307 for (
int phase = 0; phase < np; ++phase) {
308 const int pos = pu.phase_pos[phase];
309 const ADB phase_fraction = wellOps().topseg2w * (wbq[phase] / wbqt);
310 cmix_s[phase] = wellOps().w2p * aliveWells_selector.select(phase_fraction,
ADB::constant(compi.col(pos)));
315 const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
317 for (
int phase = 0; phase < np; ++phase) {
318 ADB tmp = cmix_s[phase];
319 if (phase == Oil && (*active_)[Gas]) {
320 const int gaspos = pu.phase_pos[Gas];
321 tmp = (tmp - rv_perfcells * cmix_s[gaspos]) / d;
323 if (phase == Gas && (*active_)[Oil]) {
324 const int oilpos = pu.phase_pos[Oil];
325 tmp = (tmp - rs_perfcells * cmix_s[oilpos]) / d;
327 volumeRatio += tmp / b_perfcells[phase];
331 ADB cqt_is = cqt_i/volumeRatio;
334 for (
int phase = 0; phase < np; ++phase) {
335 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
344 template <
class SolutionState,
class WellState>
347 updatePerfPhaseRatesAndPressures(
const std::vector<ADB>& cq_s,
348 const SolutionState& state,
351 if ( !localWellsActive() ) {
355 const int np = numPhases();
356 const int nw = numWells();
358 Vector cq =
superset(cq_s[0].value(), Span(nperf_total_, np, 0), nperf_total_ * np);
359 for (
int phase = 1; phase < np; ++phase) {
360 cq +=
superset(cq_s[phase].value(), Span(nperf_total_, np, phase), nperf_total_ * np);
362 xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total_ * np);
367 xw.perfPress().resize(nperf_total_, -1.e100);
369 const Vector& cdp = well_perforation_pressure_diffs_;
370 int start_segment = 0;
371 int start_perforation = 0;
372 for (
int i = 0; i < nw; ++i) {
373 WellMultiSegmentConstPtr well = wells_multisegment_[i];
374 const int nperf = well->numberOfPerforations();
375 const int nseg = well->numberOfSegments();
376 if (well->isMultiSegmented()) {
377 start_segment += nseg;
378 start_perforation += nperf;
381 const Vector cdp_well =
subset(cdp, Span(nperf, 1, start_perforation));
382 const ADB segp =
subset(state.segp, Span(nseg, 1, start_segment));
383 const Vector perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well;
384 std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]);
386 start_segment += nseg;
387 start_perforation += nperf;
389 assert(start_segment == nseg_total_);
390 assert(start_perforation == nperf_total_);
397 template <
class SolutionState>
400 computeSegmentFluidProperties(
const SolutionState& state)
402 const int np = numPhases();
403 const int nw = msWells().size();
404 const int nseg_total = nseg_total_;
406 if ( !wellOps().has_multisegment_wells ){
409 well_segment_densities_ =
ADB::constant(Vector::Zero(nseg_total));
410 segment_mass_flow_rates_ =
ADB::constant(Vector::Zero(nseg_total));
411 segment_viscosities_ =
ADB::constant(Vector::Zero(nseg_total));
412 for (
int phase = 0; phase < np; ++phase) {
413 segment_comp_surf_volume_current_[phase] =
ADB::constant(Vector::Zero(nseg_total));
414 segmentCompSurfVolumeInitial()[phase] = Vector::Zero(nseg_total);
428 std::vector<int> segment_cells;
429 segment_cells.reserve(nseg_total);
430 for (
int w = 0; w < nw; ++w) {
431 const std::vector<int>& segment_cells_well = msWells()[w]->segmentCells();
432 segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end());
434 assert(
int(segment_cells.size()) == nseg_total);
436 const ADB segment_temp =
subset(state.temperature, segment_cells);
439 const ADB& segment_press = state.segp;
442 std::vector<PhasePresence> segment_cond(nseg_total);
443 for (
int s = 0; s < nseg_total; ++s) {
444 segment_cond[s] = (*phase_condition_)[segment_cells[s]];
448 std::vector<ADB> mu_seg(np,
ADB::null());
452 if (pu.phase_used[Water]) {
453 b_seg[pu.phase_pos[Water]] = fluid_->
bWat(segment_press, segment_temp, segment_cells);
454 mu_seg[pu.phase_pos[Water]] = fluid_->
muWat(segment_press, segment_temp, segment_cells);
456 assert((*active_)[Oil]);
457 const ADB segment_so =
subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
458 if (pu.phase_used[Oil]) {
459 const ADB segment_rs =
subset(state.rs, segment_cells);
460 b_seg[pu.phase_pos[Oil]] = fluid_->
bOil(segment_press, segment_temp, segment_rs,
461 segment_cond, segment_cells);
463 rsmax_seg = fluid_->
rsSat(segment_press, segment_so, segment_cells);
464 mu_seg[pu.phase_pos[Oil]] = fluid_->
muOil(segment_press, segment_temp, segment_rs,
465 segment_cond, segment_cells);
467 assert((*active_)[Gas]);
468 if (pu.phase_used[Gas]) {
469 const ADB segment_rv =
subset(state.rv, segment_cells);
470 b_seg[pu.phase_pos[Gas]] = fluid_->
bGas(segment_press, segment_temp, segment_rv,
471 segment_cond, segment_cells);
473 rvmax_seg = fluid_->
rvSat(segment_press, segment_so, segment_cells);
474 mu_seg[pu.phase_pos[Gas]] = fluid_->
muGas(segment_press, segment_temp, segment_rv,
475 segment_cond, segment_cells);
479 ADB tot_surface_rate =
ADB::constant(Vector::Zero(nseg_total));
481 for (
int phase = 0; phase < np; ++phase) {
482 segqs[phase] =
subset(state.segqs, Span(nseg_total, 1, phase * nseg_total));
483 tot_surface_rate += segqs[phase];
487 std::vector<std::vector<double>> comp_frac(np, std::vector<double>(nseg_total, 0.0));
488 int start_segment = 0;
489 for (
int w = 0; w < nw; ++w) {
490 WellMultiSegmentConstPtr well = msWells()[w];
491 const int nseg = well->numberOfSegments();
492 const std::vector<double>& comp_frac_well = well->compFrac();
493 for (
int phase = 0; phase < np; ++phase) {
494 for (
int s = 0; s < nseg; ++s) {
495 comp_frac[phase][s + start_segment] = comp_frac_well[phase];
498 start_segment += nseg;
500 assert(start_segment == nseg_total);
505 for (
int phase = 0; phase < np; ++phase) {
508 mix[phase] =
ADB::constant(Eigen::Map<Vector>(comp_frac[phase].data(), nseg_total));
511 Selector<double> non_zero_tot_rate(tot_surface_rate.value(), Selector<double>::NotEqualZero);
512 for (
int phase = 0; phase < np; ++phase) {
513 mix[phase] = non_zero_tot_rate.select(segqs[phase] / tot_surface_rate, mix[phase]);
519 const int gaspos = pu.phase_pos[Gas];
520 const int oilpos = pu.phase_pos[Oil];
521 Selector<double> non_zero_mix_oilpos(mix[oilpos].value(), Selector<double>::GreaterZero);
522 Selector<double> non_zero_mix_gaspos(mix[gaspos].value(), Selector<double>::GreaterZero);
525 ADB big_values =
ADB::constant(Vector::Constant(nseg_total, 1.e100));
526 ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
527 ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
528 if ((*active_)[Oil]) {
529 Vector selectorUnderRsmax = Vector::Zero(nseg_total);
530 Vector selectorAboveRsmax = Vector::Zero(nseg_total);
531 for (
int s = 0; s < nseg_total; ++s) {
532 if (mix_gas_oil.value()[s] > rsmax_seg.value()[s]) {
533 selectorAboveRsmax[s] = 1.0;
535 selectorUnderRsmax[s] = 1.0;
538 rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
540 if ((*active_)[Gas]) {
541 Vector selectorUnderRvmax = Vector::Zero(nseg_total);
542 Vector selectorAboveRvmax = Vector::Zero(nseg_total);
543 for (
int s = 0; s < nseg_total; ++s) {
544 if (mix_oil_gas.value()[s] > rvmax_seg.value()[s]) {
545 selectorAboveRvmax[s] = 1.0;
547 selectorUnderRvmax[s] = 1.0;
550 rv = non_zero_mix_gaspos.select(selectorAboveRvmax * rvmax_seg + selectorUnderRvmax * mix_oil_gas, rv);
555 for (
int phase = 0; phase < np; ++phase) {
556 x[phase] = mix[phase];
558 if ((*active_)[Gas] && (*active_)[Oil]) {
559 x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
560 x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
565 for (
int phase = 0; phase < np; ++phase) {
566 volrat += x[phase] / b_seg[phase];
571 for (
int phase = 0; phase < np; ++phase) {
572 const Vector surface_density = fluid_->
surfaceDensity(phase, segment_cells);
573 dens += surface_density * mix[phase];
575 well_segment_densities_ = dens / volrat;
578 assert(np ==
int(segment_comp_surf_volume_current_.size()));
579 const ADB segment_surface_volume = segvdt_ / volrat;
580 for (
int phase = 0; phase < np; ++phase) {
581 segment_comp_surf_volume_current_[phase] = segment_surface_volume * mix[phase];
585 segment_mass_flow_rates_ =
ADB::constant(Vector::Zero(nseg_total));
586 for (
int phase = 0; phase < np; ++phase) {
588 const Vector surface_density = fluid_->
surfaceDensity(phase, segment_cells);
589 segment_mass_flow_rates_ += surface_density * segqs[phase];
593 segment_viscosities_ =
ADB::constant(Vector::Zero(nseg_total));
594 for (
int phase = 0; phase < np; ++phase) {
595 segment_viscosities_ += x[phase] * mu_seg[phase];
603 template <
class SolutionState>
606 addWellFluxEq(
const std::vector<ADB>& cq_s,
607 const SolutionState& state,
608 LinearisedBlackoilResidual& residual)
610 if ( !localWellsActive() ) {
624 const int np = numPhases();
625 const int nseg_total = nseg_total_;
627 ADB segqs = state.segqs;
629 std::vector<ADB> segment_volume_change_dt(np,
ADB::null());
630 for (
int phase = 0; phase < np; ++phase) {
631 if ( wellOps().has_multisegment_wells ) {
633 segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
634 segmentCompSurfVolumeInitial()[phase];
638 if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
639 assert(segment_volume_change_dt[phase].numBlocks() > 2);
640 assert(segqs.numBlocks() == 2);
641 segment_volume_change_dt[phase] = wellhelpers::onlyWellDerivs(segment_volume_change_dt[phase]);
642 assert(segment_volume_change_dt[phase].numBlocks() == 2);
645 const ADB cq_s_seg = wellOps().p2s * cq_s[phase];
646 const ADB segqs_phase =
subset(segqs, Span(nseg_total, 1, phase * nseg_total));
647 segqs -=
superset(cq_s_seg + wellOps().s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
648 Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
650 segqs -=
superset(wellOps().p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
654 residual.well_flux_eq = segqs;
661 template <
class SolutionState,
class WellState>
664 addWellControlEq(
const SolutionState& state,
666 const Vector& aliveWells,
667 LinearisedBlackoilResidual& residual)
672 if( msWells().empty() )
return;
674 const int np = numPhases();
675 const int nw = msWells().size();
676 const int nseg_total = nseg_total_;
682 if ((*active_)[Water]) {
683 aqua +=
subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
685 if ((*active_)[Oil]) {
686 liquid +=
subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
688 if ((*active_)[Gas]) {
689 vapour +=
subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
695 Vector rho_v = Vector::Zero(nw);
696 Vector vfp_ref_depth_v = Vector::Zero(nw);
699 Vector bhp_targets = Vector::Zero(nw);
700 Vector rate_targets = Vector::Zero(nw);
701 Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
705 std::vector<int> bhp_well_elems;
706 std::vector<int> rate_well_elems;
708 std::vector<int> bhp_top_elems;
709 std::vector<int> rate_top_elems;
710 std::vector<int> rate_top_phase_elems;
711 std::vector<int> others_elems;
715 int start_segment = 0;
716 for (
int w = 0; w < nw; ++w) {
717 const struct WellControls* wc = msWells()[w]->wellControls();
722 const int current = xw.currentControls()[w];
724 const int nseg = msWells()[w]->numberOfSegments();
726 switch (well_controls_iget_type(wc, current)) {
729 bhp_well_elems.push_back(w);
730 bhp_top_elems.push_back(start_segment);
731 bhp_targets(w) = well_controls_iget_target(wc, current);
732 rate_targets(w) = -1e100;
733 for (
int p = 0; p < np; ++p) {
734 rate_top_phase_elems.push_back(np * start_segment + p);
741 OPM_THROW(std::runtime_error,
"THP control is not implemented for multi-sgement wells yet!!");
748 rate_well_elems.push_back(w);
749 rate_top_elems.push_back(start_segment);
750 for (
int p = 0; p < np; ++p) {
751 rate_top_phase_elems.push_back(np * start_segment + p);
757 const double*
const distr =
758 well_controls_iget_distr(wc, current);
760 for (
int p = 0; p < np; ++p) {
761 rate_distr.insert(w, p*nw + w) = distr[p];
764 bhp_targets(w) = -1.0e100;
765 rate_targets(w) = well_controls_iget_target(wc, current);
771 for (
int i = 1; i < nseg; ++i) {
772 others_elems.push_back(i + start_segment);
774 start_segment += nseg;
779 const ADB bhp_residual =
subset(state.segp, bhp_top_elems) -
subset(bhp_targets, bhp_well_elems);
780 const ADB rate_residual =
subset(rate_distr *
subset(state.segqs, rate_top_phase_elems) - rate_targets, rate_well_elems);
782 ADB others_residual =
ADB::constant(Vector::Zero(nseg_total));
784 if ( wellOps().has_multisegment_wells ) {
787 ADB wspd = (state.segp.numBlocks() == 2)
788 ? wellhelpers::onlyWellDerivs(well_segment_pressures_delta_)
789 : well_segment_pressures_delta_;
791 others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp + wspd);
793 others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp);
798 ADB well_eq_topsegment =
subset(
superset(bhp_residual, bhp_top_elems, nseg_total) +
799 superset(rate_residual, rate_top_elems, nseg_total), top_well_segments_);
804 Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
805 for (
int w = 0; w < nw; ++w) {
806 for (
int phase = 0; phase < np; ++phase) {
807 rate_summer.insert(w, phase*nw + w) = 1.0;
810 Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
815 well_eq_topsegment = alive_selector.select(well_eq_topsegment, rate_summer *
subset(state.segqs, rate_top_phase_elems));
820 residual.well_eq =
superset(well_eq_topsegment, top_well_segments_, nseg_total) +
828 template <
class WellState>
831 updateWellControls(WellState& xw)
const 833 wellhelpers::WellSwitchingLogger logger;
835 if( msWells().empty() ) return ;
839 const int np = numPhases();
840 const int nw = msWells().size();
841 for (
int w = 0; w < nw; ++w) {
842 const WellControls* wc = msWells()[w]->wellControls();
846 int current = xw.currentControls()[w];
850 const int nwc = well_controls_get_num(wc);
852 for (; ctrl_index < nwc; ++ctrl_index) {
853 if (ctrl_index == current) {
859 if (wellhelpers::constraintBroken(
860 xw.bhp(), xw.thp(), xw.wellRates(),
861 w, np, msWells()[w]->wellType(), wc, ctrl_index)) {
867 if (ctrl_index != nwc) {
870 logger.wellSwitched(msWells()[w]->name(),
871 well_controls_iget_type(wc, current),
872 well_controls_iget_type(wc, ctrl_index));
873 xw.currentControls()[w] = ctrl_index;
874 current = xw.currentControls()[w];
882 const double target = well_controls_iget_target(wc, current);
883 const double* distr = well_controls_iget_distr(wc, current);
884 switch (well_controls_iget_type(wc, current)) {
886 xw.bhp()[w] = target;
887 xw.segPress()[top_well_segments_[w]] = target;
891 OPM_THROW(std::runtime_error,
"THP control is not implemented for multi-sgement wells yet!!");
902 for (
int phase = 0; phase < np; ++phase) {
903 if (distr[phase] > 0.0) {
904 xw.wellRates()[np * w + phase] = target * distr[phase];
907 xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
913 if (wellCollection()->groupControlActive()) {
915 WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
918 if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
920 well_node.setIndividualControl(
false);
923 well_node.setIndividualControl(
true);
937 template <
class SolutionState,
class WellState>
940 computeWellConnectionPressures(
const SolutionState& state,
942 const std::vector<ADB>& kr_adb,
943 const std::vector<ADB>& fluid_density)
945 if( ! wellsActive() ) return ;
951 const int nperf_total = nperf_total_;
952 const int nw = numWells();
954 const std::vector<int>& well_cells = wellOps().well_cells;
956 well_perforation_densities_ = Vector::Zero(nperf_total);
958 const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf_total);
960 Vector avg_press = perf_press * 0.0;
965 int start_segment = 0;
966 for (
int w = 0; w < nw; ++w) {
967 const int nseg = wells_multisegment_[w]->numberOfSegments();
968 if (wells_multisegment_[w]->isMultiSegmented()) {
970 start_segment += nseg;
974 std::string well_name(wells_multisegment_[w]->name());
975 typedef typename WellState::SegmentedWellMapType::const_iterator const_iterator;
976 const_iterator it_well = xw.segmentedWellMap().find(well_name);
977 assert(it_well != xw.segmentedWellMap().end());
979 const int start_perforation = (*it_well).second.start_perforation;
980 const int end_perforation = start_perforation + (*it_well).second.number_of_perforations;
981 for (
int perf = start_perforation; perf < end_perforation; ++perf) {
982 const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1];
983 const double p_avg = (perf_press[perf] + p_above)/2;
984 avg_press[perf] = p_avg;
986 start_segment += nseg;
988 assert(start_segment == xw.numSegments());
991 const ADB perf_temp =
subset(state.temperature, well_cells);
997 std::vector<PhasePresence> perf_cond(nperf_total);
998 for (
int perf = 0; perf < nperf_total; ++perf) {
999 perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
1001 const PhaseUsage& pu = fluid_->phaseUsage();
1002 DataBlock b(nperf_total, pu.num_phases);
1003 std::vector<double> rsmax_perf(nperf_total, 0.0);
1004 std::vector<double> rvmax_perf(nperf_total, 0.0);
1005 if (pu.phase_used[BlackoilPhases::Aqua]) {
1006 const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
1007 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
1009 assert((*active_)[Oil]);
1010 const Vector perf_so =
subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
1011 if (pu.phase_used[BlackoilPhases::Liquid]) {
1012 const ADB perf_rs =
subset(state.rs, well_cells);
1013 const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
1014 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
1016 rsmax_perf.assign(rssat.data(), rssat.data() + nperf_total);
1018 if (pu.phase_used[BlackoilPhases::Vapour]) {
1019 const ADB perf_rv =
subset(state.rv, well_cells);
1020 const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
1021 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
1023 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf_total);
1026 std::vector<double> b_perf(b.data(), b.data() + nperf_total * pu.num_phases);
1028 const Vector& perfcelldepth = perf_cell_depth_;
1029 std::vector<double> perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total);
1034 Vector rho =
superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases);
1035 for (
int phase = 1; phase < pu.num_phases; ++phase) {
1036 rho +=
superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases);
1038 std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases);
1041 std::vector<double> cd =
1043 wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
1044 b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1047 std::vector<double> cdp =
1049 wells(), perf_cell_depth, cd, gravity_);
1052 well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf_total);
1053 well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf_total);
1055 if ( !wellOps().has_multisegment_wells ) {
1056 well_perforation_cell_densities_ = Vector::Zero(nperf_total);
1057 well_perforation_cell_pressure_diffs_ = Vector::Zero(nperf_total);
1064 size_t temp_size = kr_adb.size();
1065 std::vector<Vector> perf_kr;
1066 for(
size_t i = 0; i < temp_size; ++i) {
1068 const Vector kr_phase = (
subset(kr_adb[i], well_cells)).value();
1069 perf_kr.push_back(kr_phase);
1076 for (
int i = 0; i < nperf_total; ++i) {
1078 int np = perf_kr.size();
1079 for (
int p = 0; p < np; ++p) {
1080 sum_kr += perf_kr[p][i];
1083 for (
int p = 0; p < np; ++p) {
1084 perf_kr[p][i] /= sum_kr;
1088 Vector rho_avg_perf = Vector::Constant(nperf_total, 0.0);
1090 for (
int phaseIdx = 0; phaseIdx < fluid_->numPhases(); ++phaseIdx) {
1093 const Vector rho_perf =
subset(fluid_density[phaseIdx], well_cells).value();
1095 rho_avg_perf += rho_perf * perf_kr[phaseIdx];
1098 well_perforation_cell_densities_ = Eigen::Map<const Vector>(rho_avg_perf.data(), nperf_total);
1100 well_perforation_cell_pressure_diffs_ = gravity_ * well_perforation_cell_densities_ * perf_cell_depth_diffs_;
1107 template <
class SolutionState>
1110 variableStateExtractWellsVars(
const std::vector<int>& indices,
1111 std::vector<ADB>& vars,
1112 SolutionState& state)
const 1118 state.segqs = std::move(vars[indices[Qs]]);
1121 state.segp = std::move(vars[indices[Bhp]]);
1126 const int np = num_phases_;
1127 const int ns = nseg_total_;
1128 const int nw = numWells();
1130 for (
int phase = 0; phase < np; ++phase) {
1132 ADB segqs_phase =
subset(state.segqs, Span(ns, 1, ns*phase));
1134 ADB wellqs_phase =
subset(segqs_phase, topWellSegments());
1136 state.qs +=
superset(wellqs_phase, Span(nw, 1, nw * phase), nw * np);
1138 state.bhp =
subset(state.segp, topWellSegments());
1145 template <
class WellState>
1148 variableWellStateInitials(
const WellState& xw,
1149 std::vector<Vector>& vars0)
const 1152 if ( wells_multisegment_.size() > 0 )
1155 const int nseg = xw.numSegments();
1156 const int np = xw.numPhases();
1159 const DataBlock segrates = Eigen::Map<const DataBlock>(& xw.segPhaseRates()[0], nseg, np).transpose();
1161 const Vector segqs = Eigen::Map<const Vector>(segrates.data(), nseg * np);
1162 vars0.push_back(segqs);
1165 const Vector segp = Eigen::Map<const Vector>(& xw.segPress()[0], xw.segPress().size());
1166 vars0.push_back(segp);
1171 vars0.push_back(Vector());
1172 vars0.push_back(Vector());
1178 #endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED ADB muOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const
Oil viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:314
ADB bOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const
Oil formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:493
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
V surfaceDensity(const int phaseIdx, const Cells &cells) const
Densities of stock components at surface conditions.
Definition: BlackoilPropsAdFromDeck.cpp:242
ADB muWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:263
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
ADB rvSat(const ADB &po, const Cells &cells) const
Condensation curve for Rv as function of oil pressure.
Definition: BlackoilPropsAdFromDeck.cpp:688
V rsSat(const V &po, const Cells &cells) const
Bubble point curve for Rs as function of oil pressure.
Definition: GridHelpers.cpp:27
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231
ADB bWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:446
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
static std::vector< double > computeConnectionPressureDelta(const Wells &wells, const std::vector< double > &z_perf, const std::vector< double > &dens_perf, const double gravity)
Compute pressure deltas.
Definition: WellDensitySegmented.cpp:271
static std::vector< double > computeConnectionDensities(const Wells &wells, const PhaseUsage &phase_usage, const std::vector< double > &perfComponentRates, const std::vector< double > &b_perf, const std::vector< double > &rsmax_perf, const std::vector< double > &rvmax_perf, const std::vector< double > &surf_dens_perf)
Compute well segment densities Notation: N = number of perforations, C = number of components...
Definition: WellDensitySegmented.cpp:32
ADB bGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const
Gas formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:566
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
ADB muGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const
Gas viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:381
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719