22#ifndef OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
23#define OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
26#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
31#include <opm/common/Exceptions.hpp>
33#include <opm/input/eclipse/Units/Units.hpp>
45#include <fmt/format.h>
50 template<
typename TypeTag>
57 const int pvtRegionIdx,
58 const int num_components,
60 const int index_of_well,
62 :
Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
73 template<
typename TypeTag>
77 const std::vector<Scalar>& depth_arg,
79 const std::vector< Scalar >& B_avg,
80 const bool changed_to_open_this_step)
82 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
83 this->StdWellEval::init(this->perf_depth_, depth_arg, Base::has_polymermw);
90 template<
typename TypeTag>
95 const std::vector<Value>& mob,
97 const std::vector<Scalar>& Tw,
100 std::vector<Value>& cq_s,
104 auto obtain = [
this](
const Eval& value)
106 if constexpr (std::is_same_v<Value, Scalar>) {
107 static_cast<void>(
this);
108 return getValue(value);
110 return this->extendEval(value);
113 auto obtainN = [](
const auto& value)
115 if constexpr (std::is_same_v<Value, Scalar>) {
116 return getValue(value);
121 auto zeroElem = [
this]()
123 if constexpr (std::is_same_v<Value, Scalar>) {
124 static_cast<void>(
this);
127 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
131 const auto& fs = intQuants.fluidState();
132 const Value pressure = obtain(this->getPerfCellPressure(fs));
133 const Value rs = obtain(fs.Rs());
134 const Value rv = obtain(fs.Rv());
135 const Value rvw = obtain(fs.Rvw());
136 const Value rsw = obtain(fs.Rsw());
138 std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
139 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
140 if (!FluidSystem::phaseIsActive(phaseIdx)) {
143 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
144 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
146 if constexpr (has_solvent) {
147 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
150 if constexpr (has_zFraction) {
151 if (this->isInjector()) {
152 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
153 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
154 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
158 Value skin_pressure = zeroElem();
160 if (this->isInjector()) {
161 const int pskin_index = Bhp + 1 + this->numLocalPerfs() + perf;
162 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
167 std::vector<Value> cmix_s(this->numComponents(), zeroElem());
168 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
169 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
192 template<
typename TypeTag>
193 template<
class Value>
197 const Value& pressure,
203 std::vector<Value>& b_perfcells_dense,
204 const std::vector<Scalar>& Tw,
207 const Value& skin_pressure,
208 const std::vector<Value>& cmix_s,
209 std::vector<Value>& cq_s,
214 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
215 Value drawdown = pressure - well_pressure;
216 if (this->isInjector()) {
217 drawdown += skin_pressure;
221 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
222 ? Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)
224 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
225 ? Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)
227 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
228 ? Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)
236 if (!allow_cf && this->isInjector()) {
241 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
242 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
243 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
246 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
247 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
249 ratioCalc.gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw,
250 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
252 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
253 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
255 ratioCalc.gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw, this->isProducer());
259 if (!allow_cf && this->isProducer()) {
264 Value total_mob_dense = mob[0];
265 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
266 total_mob_dense += mob[componentIdx];
270 Value volumeRatio = bhp * 0.0;
272 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
273 ratioCalc.disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
274 cmix_s, b_perfcells_dense, deferred_logger);
278 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
279 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
280 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
283 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
284 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
285 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
288 if constexpr (Indices::enableSolvent) {
289 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
292 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
293 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
295 ratioCalc.gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
296 cmix_s, b_perfcells_dense,
299 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
300 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
301 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
303 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
304 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
305 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
311 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
312 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
313 Value cqt_is = cqt_i / volumeRatio;
314 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
318 if (this->isProducer()) {
319 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
320 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
322 ratioCalc.gasOilPerfRateInj(cq_s, perf_rates,
323 rv, rs, pressure, rvw,
324 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
327 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
328 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
331 ratioCalc.gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
332 pressure, deferred_logger);
339 template<
typename TypeTag>
344 const Well::InjectionControls& inj_controls,
345 const Well::ProductionControls& prod_controls,
352 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
355 this->linSys_.clear();
357 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
358 prod_controls, well_state,
359 group_state, deferred_logger);
365 template<
typename TypeTag>
370 const Well::InjectionControls& inj_controls,
371 const Well::ProductionControls& prod_controls,
377 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
378 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
380 auto& ws = well_state.
well(this->index_of_well_);
381 ws.phase_mixing_rates.fill(0.0);
384 const int np = this->number_of_phases_;
386 std::vector<RateVector> connectionRates = this->connectionRates_;
388 auto& perf_data = ws.perf_data;
389 auto& perf_rates = perf_data.phase_rates;
390 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
392 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
393 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
394 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
395 calculateSinglePerf(simulator, perf, well_state, connectionRates,
396 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
399 if constexpr (has_polymer && Base::has_polymermw) {
400 if (this->isInjector()) {
401 handleInjectivityEquations(simulator, well_state, perf,
402 water_flux_s, deferred_logger);
405 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
407 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
409 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
412 assemblePerforationEq(cq_s_effective,
415 this->primary_variables_.numWellEq(),
419 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
420 auto& perf_rate_solvent = perf_data.solvent_rates;
421 perf_rate_solvent[perf] = cq_s[componentIdx].value();
423 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
427 if constexpr (has_zFraction) {
429 assembleZFracEq(cq_s_zfrac_effective,
431 this->primary_variables_.numWellEq(),
436 this->connectionRates_ = connectionRates;
441 const auto& comm = this->parallel_well_info_.communication();
442 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
446 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
449 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
452 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
453 if (FluidSystem::numActivePhases() > 1) {
455 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
456 this->F0_[componentIdx]) * volume / dt;
458 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
460 assembleSourceEq(resWell_loc,
462 this->primary_variables_.numWellEq(),
466 const auto& summaryState = simulator.vanguard().summaryState();
467 const Schedule& schedule = simulator.vanguard().schedule();
468 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
470 assembleControlEq(well_state, group_state,
471 schedule, summaryState,
472 inj_controls, prod_controls,
473 this->primary_variables_,
474 this->connections_.rho(),
476 stopped_or_zero_target,
482 this->linSys_.invert();
484 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
491 template<
typename TypeTag>
497 std::vector<RateVector>& connectionRates,
498 std::vector<EvalWell>& cq_s,
503 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
504 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
505 const int cell_idx = this->well_cells_[perf];
506 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
507 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
508 getMobility(simulator, perf, mob, deferred_logger);
511 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
512 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
513 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
514 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
515 cq_s, perf_rates, deferred_logger);
517 auto& ws = well_state.
well(this->index_of_well_);
518 auto& perf_data = ws.perf_data;
519 if constexpr (has_polymer && Base::has_polymermw) {
520 if (this->isInjector()) {
523 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
524 water_flux_s = cq_s[water_comp_idx];
527 handleInjectivityRate(simulator, perf, cq_s);
532 if (this->isProducer()) {
533 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.
dis_gas;
534 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.
dis_gas_in_water;
535 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.
vap_oil;
536 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.
vap_wat;
537 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.
dis_gas;
538 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.
dis_gas_in_water;
539 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.
vap_oil;
540 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.
vap_wat;
543 if constexpr (has_energy) {
544 connectionRates[perf][Indices::contiEnergyEqIdx] =
545 connectionRateEnergy(cq_s, intQuants, deferred_logger);
548 if constexpr (has_polymer) {
549 std::variant<Scalar,EvalWell> polymerConcentration;
550 if (this->isInjector()) {
551 polymerConcentration = this->wpolymer();
553 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
554 intQuants.polymerViscosityCorrection());
557 [[maybe_unused]]
EvalWell cq_s_poly;
558 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
560 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
561 cq_s, polymerConcentration);
563 if constexpr (Base::has_polymermw) {
564 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
565 perf, connectionRates, deferred_logger);
569 if constexpr (has_foam) {
570 std::variant<Scalar,EvalWell> foamConcentration;
571 if (this->isInjector()) {
572 foamConcentration = this->wfoam();
574 foamConcentration = this->extendEval(intQuants.foamConcentration());
576 connectionRates[perf][Indices::contiFoamEqIdx] =
577 this->connections_.connectionRateFoam(cq_s, foamConcentration,
578 FoamModule::transportPhase(),
582 if constexpr (has_zFraction) {
583 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
584 if (this->isInjector()) {
585 solventConcentration = this->wsolvent();
587 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
588 this->extendEval(intQuants.yVolume())};
590 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
591 cq_s_zfrac_effective) =
592 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
594 solventConcentration);
597 if constexpr (has_brine) {
598 std::variant<Scalar,EvalWell> saltConcentration;
599 if (this->isInjector()) {
600 saltConcentration = this->wsalt();
602 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
605 connectionRates[perf][Indices::contiBrineEqIdx] =
606 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
611 if constexpr (has_micp) {
612 std::variant<Scalar,EvalWell> microbialConcentration;
613 std::variant<Scalar,EvalWell> oxygenConcentration;
614 std::variant<Scalar,EvalWell> ureaConcentration;
615 if (this->isInjector()) {
616 microbialConcentration = this->wmicrobes();
617 oxygenConcentration = this->woxygen();
618 ureaConcentration = this->wurea();
620 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
621 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
622 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
624 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
625 connectionRates[perf][Indices::contiOxygenEqIdx],
626 connectionRates[perf][Indices::contiUreaEqIdx]) =
627 this->connections_.connectionRatesMICP(perf_data.microbial_rates[perf],
628 perf_data.oxygen_rates[perf],
629 perf_data.urea_rates[perf],
631 microbialConcentration,
637 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
641 if (pu.has_co2_or_h2store) {
642 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
643 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
644 perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
648 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
649 const unsigned wat_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
650 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::waterPhaseIdx, Base::pvtRegionIdx() );
651 perf_data.wat_mass_rates[perf] = cq_s[wat_comp_idx].value() * rho;
657 template<
typename TypeTag>
658 template<
class Value>
663 std::vector<Value>& mob,
666 auto obtain = [
this](
const Eval& value)
668 if constexpr (std::is_same_v<Value, Scalar>) {
669 static_cast<void>(
this);
670 return getValue(value);
672 return this->extendEval(value);
676 obtain, deferred_logger);
679 if constexpr (has_polymer) {
680 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
681 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
686 if constexpr (!Base::has_polymermw) {
687 if constexpr (std::is_same_v<Value, Scalar>) {
688 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
689 for (std::size_t i = 0; i < mob.size(); ++i) {
690 mob_eval[i].setValue(mob[i]);
692 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
693 for (std::size_t i = 0; i < mob.size(); ++i) {
694 mob[i] = getValue(mob_eval[i]);
697 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
704 const Scalar bhp = this->primary_variables_.value(Bhp);
705 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
706 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
707 for (std::size_t i = 0; i < mob.size(); ++i) {
708 mob[i] *= multiplier;
714 template<
typename TypeTag>
722 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
724 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
725 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
727 const auto& summary_state = simulator.vanguard().summaryState();
728 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
729 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.
well(this->index_of_well_));
736 template<
typename TypeTag>
740 const bool stop_or_zero_rate_target,
743 const Scalar dFLimit = this->param_.dwell_fraction_max_;
744 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
745 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
748 if constexpr (Base::has_polymermw) {
749 this->primary_variables_.updateNewtonPolyMW(dwells);
752 this->primary_variables_.checkFinite(deferred_logger);
759 template<
typename TypeTag>
763 const SummaryState& summary_state,
766 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
769 if constexpr (Base::has_polymermw) {
770 this->primary_variables_.copyToWellStatePolyMW(well_state);
778 template<
typename TypeTag>
786 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
787 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
789 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
790 std::vector<Scalar> mob(this->num_components_, 0.0);
791 getMobility(simulator, perf, mob, deferred_logger);
793 const int cell_idx = this->well_cells_[perf];
794 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
795 const auto& fs = int_quantities.fluidState();
797 Scalar p_r = this->getPerfCellPressure(fs).value();
800 std::vector<Scalar> b_perf(this->num_components_);
801 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
802 if (!FluidSystem::phaseIsActive(phase)) {
805 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
806 b_perf[comp_idx] = fs.invB(phase).value();
808 if constexpr (has_solvent) {
809 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
813 const Scalar h_perf = this->connections_.pressure_diff(perf);
814 const Scalar pressure_diff = p_r - h_perf;
819 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
820 deferred_logger.
debug(
"CROSSFLOW_IPR",
821 "cross flow found when updateIPR for well " + name()
822 +
" . The connection is ignored in IPR calculations");
828 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
829 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
830 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
834 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
835 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
836 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
837 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
838 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
839 ipr_b_perf[comp_idx] += tw_mob;
843 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
844 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
845 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
846 const Scalar rs = (fs.Rs()).value();
847 const Scalar rv = (fs.Rv()).value();
849 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
850 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
852 ipr_a_perf[gas_comp_idx] += dis_gas_a;
853 ipr_a_perf[oil_comp_idx] += vap_oil_a;
855 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
856 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
858 ipr_b_perf[gas_comp_idx] += dis_gas_b;
859 ipr_b_perf[oil_comp_idx] += vap_oil_b;
862 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
863 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
864 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
867 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
868 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
871 template<
typename TypeTag>
885 auto rates = well_state.
well(this->index_of_well_).surface_rates;
887 for (std::size_t p = 0; p < rates.size(); ++p) {
888 zero_rates &= rates[p] == 0.0;
890 auto& ws = well_state.
well(this->index_of_well_);
892 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
893 deferred_logger.
debug(msg);
905 const auto& group_state = simulator.problem().wellModel().groupState();
907 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
908 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
910 auto inj_controls = Well::InjectionControls(0);
911 auto prod_controls = Well::ProductionControls(0);
912 prod_controls.addControl(Well::ProducerCMode::BHP);
913 prod_controls.bhp_limit = well_state.
well(this->index_of_well_).bhp;
916 const auto cmode = ws.production_cmode;
917 ws.production_cmode = Well::ProducerCMode::BHP;
918 const double dt = simulator.timeStepSize();
919 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
921 const size_t nEq = this->primary_variables_.numWellEq();
925 for (
size_t i=0; i < nEq; ++i){
931 x_well[0].resize(nEq);
932 this->linSys_.solve(rhs, x_well);
934 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
935 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
936 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
937 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
939 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
941 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
944 ws.production_cmode = cmode;
947 template<
typename TypeTag>
954 const auto& summaryState = simulator.vanguard().summaryState();
958 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
959 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
962 Scalar total_ipr_mass_rate = 0.0;
963 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
965 if (!FluidSystem::phaseIsActive(phaseIdx)) {
969 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
970 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
972 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
973 total_ipr_mass_rate += ipr_rate * rho;
975 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
976 this->operability_status_.operable_under_only_bhp_limit =
false;
980 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
984 std::vector<Scalar> well_rates_bhp_limit;
985 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
987 this->adaptRatesForVFP(well_rates_bhp_limit);
988 const Scalar thp_limit = this->getTHPConstraint(summaryState);
991 this->connections_.rho(),
992 this->getALQ(well_state),
995 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
996 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1007 this->operability_status_.operable_under_only_bhp_limit =
true;
1008 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1016 template<
typename TypeTag>
1023 const auto& summaryState = simulator.vanguard().summaryState();
1024 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1025 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1028 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1031 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1032 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1034 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1035 if (this->isProducer() && *obtain_bhp < thp_limit) {
1036 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1037 +
" bars is SMALLER than thp limit "
1039 +
" bars as a producer for well " + name();
1040 deferred_logger.
debug(msg);
1042 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1043 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1044 +
" bars is LARGER than thp limit "
1046 +
" bars as a injector for well " + name();
1047 deferred_logger.
debug(msg);
1050 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1051 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1052 if (!this->wellIsStopped()) {
1053 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1054 deferred_logger.
debug(
" could not find bhp value at thp limit "
1056 +
" bar for well " + name() +
", the well might need to be closed ");
1065 template<
typename TypeTag>
1070 bool all_drawdown_wrong_direction =
true;
1072 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1073 const int cell_idx = this->well_cells_[perf];
1074 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1075 const auto& fs = intQuants.fluidState();
1077 const Scalar pressure = this->getPerfCellPressure(fs).value();
1078 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1081 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1082 const Scalar drawdown = pressure - well_pressure;
1087 if ( (drawdown < 0. && this->isInjector()) ||
1088 (drawdown > 0. && this->isProducer()) ) {
1089 all_drawdown_wrong_direction =
false;
1094 const auto& comm = this->parallel_well_info_.communication();
1095 if (comm.size() > 1)
1097 all_drawdown_wrong_direction =
1098 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1101 return all_drawdown_wrong_direction;
1107 template<
typename TypeTag>
1112 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1118 template<
typename TypeTag>
1124 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1126 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1128 return model.intensiveQuantities(cell_idx, 0)
1129 .fluidState().temperature(phase_idx).value();
1133 [&model = simulator.model()](
int cell_idx)
1135 return model.intensiveQuantities(cell_idx, 0)
1136 .fluidState().saltConcentration().value();
1140 [&model = simulator.model()](
int cell_idx)
1142 return model.intensiveQuantities(cell_idx, 0)
1143 .fluidState().pvtRegionIndex();
1147 if constexpr (Indices::enableSolvent) {
1148 prop_func.solventInverseFormationVolumeFactor =
1149 [&model = simulator.model()](
int cell_idx)
1151 return model.intensiveQuantities(cell_idx, 0)
1152 .solventInverseFormationVolumeFactor().value();
1155 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1157 return model.intensiveQuantities(cell_idx, 0)
1158 .solventRefDensity();
1162 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1169 template<
typename TypeTag>
1174 const std::vector<Scalar>& B_avg,
1176 const bool relax_tolerance)
const
1180 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1182 Scalar tol_wells = this->param_.tolerance_wells_;
1184 constexpr Scalar stopped_factor = 1.e-4;
1186 constexpr Scalar dynamic_thp_factor = 1.e-1;
1187 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1188 tol_wells = tol_wells*stopped_factor;
1189 }
else if (this->getDynamicThpLimit()) {
1190 tol_wells = tol_wells*dynamic_thp_factor;
1193 std::vector<Scalar> res;
1196 this->param_.max_residual_allowed_,
1198 this->param_.relaxed_tolerance_flow_well_,
1200 this->wellIsStopped(),
1204 checkConvergenceExtraEqs(res, report);
1213 template<
typename TypeTag>
1221 auto fluidState = [&simulator,
this](
const int perf)
1223 const auto cell_idx = this->well_cells_[perf];
1224 return simulator.model()
1225 .intensiveQuantities(cell_idx, 0).fluidState();
1228 const int np = this->number_of_phases_;
1229 auto setToZero = [np](
Scalar* x) ->
void
1231 std::fill_n(x, np, 0.0);
1234 auto addVector = [np](
const Scalar* src,
Scalar* dest) ->
void
1236 std::transform(src, src + np, dest, dest, std::plus<>{});
1239 auto& ws = well_state.
well(this->index_of_well_);
1240 auto& perf_data = ws.perf_data;
1241 auto* wellPI = ws.productivity_index.data();
1242 auto* connPI = perf_data.prod_index.data();
1246 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1247 auto subsetPerfID = 0;
1249 for (
const auto& perf : *this->perf_data_) {
1250 auto allPerfID = perf.ecl_index;
1252 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) ->
Scalar
1257 std::vector<Scalar> mob(this->num_components_, 0.0);
1258 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1260 const auto& fs = fluidState(subsetPerfID);
1263 if (this->isInjector()) {
1264 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1265 mob, connPI, deferred_logger);
1268 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1271 addVector(connPI, wellPI);
1278 const auto& comm = this->parallel_well_info_.communication();
1279 if (comm.size() > 1) {
1280 comm.sum(wellPI, np);
1283 assert ((
static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1284 "Internal logic error in processing connections for PI/II");
1289 template<
typename TypeTag>
1299 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1304 [&model = simulator.model()](
const int cell,
1305 const std::vector<int>& phases,
1306 std::vector<Scalar>& mob)
1308 const auto& iq = model.intensiveQuantities(cell, 0);
1310 std::transform(phases.begin(), phases.end(), mob.begin(),
1311 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1316 [&model = simulator.model()](
const int cell,
1317 const std::vector<int>& phases,
1318 std::vector<Scalar>& rho)
1320 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1322 std::transform(phases.begin(), phases.end(), rho.begin(),
1323 [&fs](
const int phase) { return fs.density(phase).value(); });
1327 const auto stopped_or_zero_rate_target = this->
1328 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1331 .computeProperties(stopped_or_zero_rate_target, well_state,
1332 prop_func, props, deferred_logger);
1339 template<
typename TypeTag>
1346 const auto props = computePropertiesForWellConnectionPressures
1347 (simulator, well_state);
1349 computeWellConnectionDensitesPressures(simulator, well_state,
1350 props, deferred_logger);
1357 template<
typename TypeTag>
1364 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1369 dx_well[0].resize(this->primary_variables_.numWellEq());
1370 this->linSys_.solve( dx_well);
1372 updateWellState(simulator, dx_well, well_state, deferred_logger);
1379 template<
typename TypeTag>
1386 updatePrimaryVariables(simulator, well_state, deferred_logger);
1387 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1388 this->computeAccumWell();
1393 template<
typename TypeTag>
1398 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1400 if (this->param_.matrix_add_well_contributions_)
1406 this->linSys_.apply(x, Ax);
1412 template<
typename TypeTag>
1417 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1419 this->linSys_.apply(r);
1425 template<
typename TypeTag>
1433 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1436 xw[0].resize(this->primary_variables_.numWellEq());
1438 this->linSys_.recoverSolutionWell(x, xw);
1439 updateWellState(simulator, xw, well_state, deferred_logger);
1445 template<
typename TypeTag>
1450 std::vector<Scalar>& well_flux,
1454 const int np = this->number_of_phases_;
1455 well_flux.resize(np, 0.0);
1457 const bool allow_cf = this->getAllowCrossFlow();
1459 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1460 const int cell_idx = this->well_cells_[perf];
1461 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1463 std::vector<Scalar> mob(this->num_components_, 0.);
1464 getMobility(simulator, perf, mob, deferred_logger);
1465 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1466 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1467 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1469 std::vector<Scalar> cq_s(this->num_components_, 0.);
1471 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1472 cq_s, perf_rates, deferred_logger);
1474 for(
int p = 0; p < np; ++p) {
1475 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1479 if constexpr (has_solvent) {
1481 assert(pu.phase_used[Gas]);
1482 const int gas_pos = pu.phase_pos[Gas];
1483 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1486 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1491 template<
typename TypeTag>
1496 std::vector<Scalar>& well_flux,
1507 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1508 const auto& group_state = simulator.problem().wellModel().groupState();
1511 const auto& summary_state = simulator.vanguard().summaryState();
1512 auto inj_controls = well_copy.
well_ecl_.isInjector()
1513 ? well_copy.
well_ecl_.injectionControls(summary_state)
1514 : Well::InjectionControls(0);
1515 auto prod_controls = well_copy.
well_ecl_.isProducer()
1516 ? well_copy.
well_ecl_.productionControls(summary_state) :
1517 Well::ProductionControls(0);
1520 auto& ws = well_state_copy.
well(this->index_of_well_);
1522 inj_controls.bhp_limit = bhp;
1523 ws.injection_cmode = Well::InjectorCMode::BHP;
1525 prod_controls.bhp_limit = bhp;
1526 ws.production_cmode = Well::ProducerCMode::BHP;
1531 const int np = this->number_of_phases_;
1532 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1533 for (
int phase = 0; phase < np; ++phase){
1534 well_state_copy.
wellRates(this->index_of_well_)[phase]
1535 = sign * ws.well_potentials[phase];
1540 const double dt = simulator.timeStepSize();
1541 const bool converged = well_copy.
iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1543 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1544 " potentials are computed based on unconverged solution";
1545 deferred_logger.
debug(msg);
1555 template<
typename TypeTag>
1556 std::vector<typename StandardWell<TypeTag>::Scalar>
1562 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1563 const auto& summary_state = simulator.vanguard().summaryState();
1565 const auto& well = this->well_ecl_;
1566 if (well.isInjector()){
1567 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1568 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1569 if (bhp_at_thp_limit) {
1570 const Scalar bhp = std::min(*bhp_at_thp_limit,
1571 static_cast<Scalar>(controls.bhp_limit));
1572 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1574 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1575 "Failed in getting converged thp based potential calculation for well "
1576 + name() +
". Instead the bhp based value is used");
1577 const Scalar bhp = controls.bhp_limit;
1578 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1581 computeWellRatesWithThpAlqProd(
1582 simulator, summary_state,
1583 deferred_logger, potentials, this->getALQ(well_state)
1590 template<
typename TypeTag>
1595 std::vector<Scalar>& well_potentials,
1605 const auto& group_state = simulator.problem().wellModel().groupState();
1606 auto& ws = well_state_copy.
well(this->index_of_well_);
1609 const auto& summary_state = simulator.vanguard().summaryState();
1610 auto inj_controls = well_copy.
well_ecl_.isInjector()
1611 ? well_copy.
well_ecl_.injectionControls(summary_state)
1612 : Well::InjectionControls(0);
1613 auto prod_controls = well_copy.
well_ecl_.isProducer()
1614 ? well_copy.
well_ecl_.productionControls(summary_state) :
1615 Well::ProductionControls(0);
1621 const int num_perf = ws.perf_data.size();
1622 for (
int perf = 0; perf < num_perf; ++perf) {
1626 const int np = this->number_of_phases_;
1627 bool trivial =
true;
1628 for (
int phase = 0; phase < np; ++phase){
1629 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1633 for (
int phase = 0; phase < np; ++phase) {
1634 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1639 const double dt = simulator.timeStepSize();
1641 bool converged =
false;
1642 if (this->well_ecl_.isProducer()) {
1643 converged = well_copy.
solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1645 converged = well_copy.
iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1649 well_potentials.clear();
1650 well_potentials.resize(np, 0.0);
1651 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1652 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1654 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1658 if constexpr (has_solvent) {
1660 assert(pu.phase_used[Gas]);
1661 const int gas_pos = pu.phase_pos[Gas];
1663 well_potentials[gas_pos] += rate.value();
1669 template<
typename TypeTag>
1673 const SummaryState &summary_state,
1675 std::vector<Scalar>& potentials,
1679 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1680 simulator, summary_state, alq, deferred_logger,
true);
1681 if (bhp_at_thp_limit) {
1682 const auto& controls = this->well_ecl_.productionControls(summary_state);
1683 bhp = std::max(*bhp_at_thp_limit,
1684 static_cast<Scalar>(controls.bhp_limit));
1685 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1688 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1689 "Failed in getting converged thp based potential calculation for well "
1690 + name() +
". Instead the bhp based value is used");
1691 const auto& controls = this->well_ecl_.productionControls(summary_state);
1692 bhp = controls.bhp_limit;
1693 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1698 template<
typename TypeTag>
1702 const SummaryState& summary_state,
1704 std::vector<Scalar>& potentials,
1708 computeWellRatesAndBhpWithThpAlqProd(simulator,
1715 template<
typename TypeTag>
1720 std::vector<Scalar>& well_potentials,
1723 const auto [compute_potential, bhp_controlled_well] =
1726 if (!compute_potential) {
1730 bool converged_implicit =
false;
1734 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1735 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1737 if (!converged_implicit) {
1739 const auto& summaryState = simulator.vanguard().summaryState();
1740 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1750 const auto& ws = well_state.
well(this->index_of_well_);
1751 if (this->isInjector())
1752 bhp = std::max(ws.bhp, bhp);
1754 bhp = std::min(ws.bhp, bhp);
1756 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1757 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1760 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1764 this->checkNegativeWellPotentials(well_potentials,
1765 this->param_.check_well_operability_,
1775 template<
typename TypeTag>
1779 const int openConnIdx)
const
1781 return (openConnIdx < 0)
1783 : this->connections_.rho(openConnIdx);
1790 template<
typename TypeTag>
1797 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1799 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1800 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1803 if constexpr (Base::has_polymermw) {
1804 this->primary_variables_.updatePolyMW(well_state);
1807 this->primary_variables_.checkFinite(deferred_logger);
1813 template<
typename TypeTag>
1818 return this->connections_.rho();
1824 template<
typename TypeTag>
1829 std::vector<EvalWell>& mob,
1832 const int cell_idx = this->well_cells_[perf];
1833 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1834 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1838 if (this->isInjector()) {
1840 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1841 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1842 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1845 if (PolymerModule::hasPlyshlog()) {
1848 if (this->isInjector() && this->wpolymer() == 0.) {
1853 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1854 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1856 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1858 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1859 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1860 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1861 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1862 perf_rates, deferred_logger);
1864 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1865 const auto& material_law_manager = simulator.problem().materialLawManager();
1866 const auto& scaled_drainage_info =
1867 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1868 const Scalar swcr = scaled_drainage_info.Swcr;
1869 const EvalWell poro = this->extendEval(int_quant.porosity());
1870 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1872 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1873 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1874 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1876 if (PolymerModule::hasShrate()) {
1879 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1881 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1882 int_quant.pvtRegionIndex(),
1885 mob[waterCompIdx] /= shear_factor;
1889 template<
typename TypeTag>
1893 this->linSys_.extract(jacobian);
1897 template <
typename TypeTag>
1901 const int pressureVarIndex,
1902 const bool use_well_weights,
1905 this->linSys_.extractCPRPressureMatrix(jacobian,
1916 template<
typename TypeTag>
1923 if constexpr (Base::has_polymermw) {
1924 const int water_table_id = this->polymerWaterTable_();
1925 if (water_table_id <= 0) {
1927 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1930 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1931 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1933 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1934 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1938 fmt::format(
"Polymermw is not activated, while injecting "
1939 "skin pressure is requested for well {}", name()),
1948 template<
typename TypeTag>
1956 if constexpr (Base::has_polymermw) {
1957 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1958 const EvalWell water_velocity_abs = abs(water_velocity);
1959 if (poly_inj_conc == 0.) {
1960 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1962 const int polymer_table_id = this->polymerTable_();
1963 if (polymer_table_id <= 0) {
1965 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1968 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1969 const Scalar reference_concentration = skprpolytable.refConcentration;
1970 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1972 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1973 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1974 if (poly_inj_conc == reference_concentration) {
1975 return sign * pskin_poly;
1978 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1979 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1980 return sign * pskin;
1983 fmt::format(
"Polymermw is not activated, while injecting "
1984 "skin pressure is requested for well {}", name()),
1993 template<
typename TypeTag>
2000 if constexpr (Base::has_polymermw) {
2001 const int table_id = this->polymerInjTable_();
2002 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2003 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2004 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2005 if (this->wpolymer() == 0.) {
2006 return molecular_weight;
2008 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2009 return molecular_weight;
2012 fmt::format(
"Polymermw is not activated, while injecting "
2013 "polymer molecular weight is requested for well {}", name()),
2022 template<
typename TypeTag>
2028 if constexpr (Base::has_polymermw) {
2029 if (!this->isInjector()) {
2033 auto& perf_water_throughput = well_state.
well(this->index_of_well_)
2034 .perf_data.water_throughput;
2036 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2037 const Scalar perf_water_vel =
2038 this->primary_variables_.value(Bhp + 1 + perf);
2042 if (perf_water_vel >
Scalar{0}) {
2043 perf_water_throughput[perf] += perf_water_vel * dt;
2053 template<
typename TypeTag>
2058 std::vector<EvalWell>& cq_s)
const
2060 const int cell_idx = this->well_cells_[perf];
2061 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2062 const auto& fs = int_quants.fluidState();
2063 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2064 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2065 const int wat_vel_index = Bhp + 1 + perf;
2066 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2070 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2076 template<
typename TypeTag>
2085 const int cell_idx = this->well_cells_[perf];
2086 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2087 const auto& fs = int_quants.fluidState();
2088 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2089 const EvalWell water_flux_r = water_flux_s / b_w;
2090 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2091 const EvalWell water_velocity = water_flux_r / area;
2092 const int wat_vel_index = Bhp + 1 + perf;
2095 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2097 const auto& ws = well_state.
well(this->index_of_well_);
2098 const auto& perf_data = ws.perf_data;
2099 const auto& perf_water_throughput = perf_data.water_throughput;
2100 const Scalar throughput = perf_water_throughput[perf];
2101 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2103 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2104 poly_conc.setValue(this->wpolymer());
2107 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2108 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2111 assembleInjectivityEq(eq_pskin,
2116 this->primary_variables_.numWellEq(),
2124 template<
typename TypeTag>
2133 if constexpr (Base::has_polymermw) {
2135 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2143 template<
typename TypeTag>
2150 std::vector<RateVector>& connectionRates,
2155 if (this->isInjector()) {
2156 const int wat_vel_index = Bhp + 1 + perf;
2157 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2158 if (water_velocity > 0.) {
2159 const auto& ws = well_state.
well(this->index_of_well_);
2160 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2161 const Scalar throughput = perf_water_throughput[perf];
2162 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2163 cq_s_polymw *= molecular_weight;
2169 }
else if (this->isProducer()) {
2170 if (cq_s_polymw < 0.) {
2171 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2178 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2186 template<
typename TypeTag>
2187 std::optional<typename StandardWell<TypeTag>::Scalar>
2191 const SummaryState& summary_state,
2194 return computeBhpAtThpLimitProdWithAlq(simulator,
2196 this->getALQ(well_state),
2201 template<
typename TypeTag>
2202 std::optional<typename StandardWell<TypeTag>::Scalar>
2205 const SummaryState& summary_state,
2208 bool iterate_if_no_solution)
const
2212 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2218 std::vector<Scalar> rates(3);
2219 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2220 this->adaptRatesForVFP(rates);
2224 Scalar max_pressure = 0.0;
2225 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2226 const int cell_idx = this->well_cells_[perf];
2227 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2228 const auto& fs = int_quants.fluidState();
2229 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2230 max_pressure = std::max(max_pressure, pressure_cell);
2235 this->connections_.rho(),
2237 this->getTHPConstraint(summary_state),
2241 auto v = frates(*bhpAtLimit);
2242 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2247 if (!iterate_if_no_solution)
2248 return std::nullopt;
2250 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2254 std::vector<Scalar> rates(3);
2255 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2256 this->adaptRatesForVFP(rates);
2263 this->connections_.rho(),
2265 this->getTHPConstraint(summary_state),
2271 auto v = frates(*bhpAtLimit);
2272 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2278 return std::nullopt;
2283 template<
typename TypeTag>
2284 std::optional<typename StandardWell<TypeTag>::Scalar>
2287 const SummaryState& summary_state,
2291 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2297 std::vector<Scalar> rates(3);
2298 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2304 this->connections_.rho(),
2315 template<
typename TypeTag>
2320 const Well::InjectionControls& inj_controls,
2321 const Well::ProductionControls& prod_controls,
2326 updatePrimaryVariables(simulator, well_state, deferred_logger);
2328 const int max_iter = this->param_.max_inner_iter_wells_;
2331 bool relax_convergence =
false;
2332 this->regularize_ =
false;
2334 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2336 if (it > this->param_.strict_inner_iter_wells_) {
2337 relax_convergence =
true;
2338 this->regularize_ =
true;
2341 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2343 converged = report.converged();
2349 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2356 }
while (it < max_iter);
2362 template<
typename TypeTag>
2367 const Well::InjectionControls& inj_controls,
2368 const Well::ProductionControls& prod_controls,
2372 const bool fixed_control ,
2373 const bool fixed_status )
2375 updatePrimaryVariables(simulator, well_state, deferred_logger);
2377 const int max_iter = this->param_.max_inner_iter_wells_;
2379 bool converged =
false;
2380 bool relax_convergence =
false;
2381 this->regularize_ =
false;
2382 const auto& summary_state = simulator.vanguard().summaryState();
2387 constexpr int min_its_after_switch = 4;
2388 int its_since_last_switch = min_its_after_switch;
2389 int switch_count= 0;
2391 const auto well_status_orig = this->wellStatus_;
2392 const auto operability_orig = this->operability_status_;
2393 auto well_status_cur = well_status_orig;
2394 int status_switch_count = 0;
2396 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2397 well_state.
well(this->index_of_well_).status == WellStatus::OPEN;
2399 const bool allow_switching =
2400 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2401 (!fixed_control || !fixed_status) && allow_open;
2403 bool changed =
false;
2404 bool final_check =
false;
2406 this->operability_status_.resetOperability();
2407 this->operability_status_.solvable =
true;
2409 its_since_last_switch++;
2410 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2411 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2412 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2413 inj_controls, prod_controls, wqTotal,
2414 deferred_logger, fixed_control, fixed_status);
2416 its_since_last_switch = 0;
2418 if (well_status_cur != this->wellStatus_) {
2419 well_status_cur = this->wellStatus_;
2420 status_switch_count++;
2423 if (!changed && final_check) {
2426 final_check =
false;
2430 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2432 if (it > this->param_.strict_inner_iter_wells_) {
2433 relax_convergence =
true;
2434 this->regularize_ =
true;
2437 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2439 converged = report.converged();
2443 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2445 its_since_last_switch = min_its_after_switch;
2452 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2454 }
while (it < max_iter);
2457 if (allow_switching){
2459 const bool is_stopped = this->wellIsStopped();
2460 if (this->wellHasTHPConstraints(summary_state)){
2461 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2462 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2464 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2468 this->wellStatus_ = well_status_orig;
2469 this->operability_status_ = operability_orig;
2470 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2471 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2472 deferred_logger.
debug(message);
2478 template<
typename TypeTag>
2479 std::vector<typename StandardWell<TypeTag>::Scalar>
2485 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2486 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2487 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2488 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2489 const int cell_idx = this->well_cells_[perf];
2490 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2491 std::vector<Scalar> mob(this->num_components_, 0.);
2492 getMobility(simulator, perf, mob, deferred_logger);
2493 std::vector<Scalar> cq_s(this->num_components_, 0.);
2494 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2495 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2496 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2498 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2499 cq_s, perf_rates, deferred_logger);
2500 for (
int comp = 0; comp < this->num_components_; ++comp) {
2501 well_q_s[comp] += cq_s[comp];
2504 const auto& comm = this->parallel_well_info_.communication();
2505 if (comm.size() > 1)
2507 comm.sum(well_q_s.data(), well_q_s.size());
2514 template <
typename TypeTag>
2515 std::vector<typename StandardWell<TypeTag>::Scalar>
2519 const int num_pri_vars = this->primary_variables_.numWellEq();
2520 std::vector<Scalar> retval(num_pri_vars);
2521 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2522 retval[ii] = this->primary_variables_.value(ii);
2531 template <
typename TypeTag>
2536 const int num_pri_vars = this->primary_variables_.numWellEq();
2537 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2538 this->primary_variables_.setValue(ii, it[ii]);
2540 return num_pri_vars;
2544 template <
typename TypeTag>
2548 const IntensiveQuantities& intQuants,
2551 auto fs = intQuants.fluidState();
2553 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2554 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2559 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2560 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2561 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2562 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2563 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2566 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2567 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2572 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2574 deferred_logger.
debug(
2575 fmt::format(
"Problematic d value {} obtained for well {}"
2576 " during calculateSinglePerf with rs {}"
2577 ", rv {}. Continue as if no dissolution (rs = 0) and"
2578 " vaporization (rv = 0) for this connection.",
2579 d, this->name(), fs.Rs(), fs.Rv()));
2580 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2582 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2583 cq_r_thermal = (cq_s[gasCompIdx] -
2584 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2585 (d * this->extendEval(fs.invB(phaseIdx)) );
2586 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2588 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2590 (d * this->extendEval(fs.invB(phaseIdx)) );
2596 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2598 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2599 fs.setTemperature(this->well_ecl_.inj_temperature());
2600 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2601 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2602 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2603 paramCache.setRegionIndex(pvtRegionIdx);
2604 paramCache.updatePhase(fs, phaseIdx);
2606 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2607 fs.setDensity(phaseIdx, rho);
2608 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2609 fs.setEnthalpy(phaseIdx, h);
2610 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2611 result += getValue(cq_r_thermal);
2612 }
else if (cq_r_thermal > 0.0) {
2613 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2614 result += Base::restrictEval(cq_r_thermal);
2617 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2618 result += Base::restrictEval(cq_r_thermal);
2622 return result * this->well_efficiency_factor_;
#define OPM_DEFLOG_THROW(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:45
#define OPM_DEFLOG_PROBLEM(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:61
Definition: ConvergenceReport.hpp:38
Definition: DeferredLogger.hpp:57
void warning(const std::string &tag, const std::string &message)
void debug(const std::string &tag, const std::string &message)
Definition: GroupState.hpp:43
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:195
Definition: RatioCalculator.hpp:38
Class handling assemble of the equation system for StandardWell.
Definition: StandardWellAssemble.hpp:43
Scalar pressure_diff(const unsigned perf) const
Returns pressure drop for a given perforation.
Definition: StandardWellConnections.hpp:105
StdWellConnections connections_
Connection level values.
Definition: StandardWellEval.hpp:104
PrimaryVariables primary_variables_
Primary variables for well.
Definition: StandardWellEval.hpp:98
Definition: StandardWell.hpp:60
virtual void init(const PhaseUsage *phase_usage_arg, const std::vector< Scalar > &depth_arg, const Scalar gravity_arg, const std::vector< Scalar > &B_avg, const bool changed_to_open_this_step) override
Definition: StandardWell_impl.hpp:76
void computeWellConnectionPressures(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1342
EvalWell wpolymermw(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1996
void computeWellPotentials(const Simulator &simulator, const WellState< Scalar > &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) override
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1718
typename StdWellEval::EvalWell EvalWell
Definition: StandardWell.hpp:119
typename StdWellEval::BVectorWell BVectorWell
Definition: StandardWell.hpp:120
void checkOperabilityUnderBHPLimit(const WellState< Scalar > &well_state, const Simulator &simulator, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:950
bool iterateWellEqWithControl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:2318
void computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1701
void addWellContributions(SparseMatrixAdapter &mat) const override
Definition: StandardWell_impl.hpp:1891
std::vector< Scalar > getPrimaryVars() const override
Definition: StandardWell_impl.hpp:2517
void updateIPRImplicit(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:874
void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator< Scalar > &wellPICalc, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1216
void updateWaterMobilityWithPolymer(const Simulator &simulator, const int perf, std::vector< EvalWell > &mob_water, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1827
bool iterateWellEqWithSwitching(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger, const bool fixed_control=false, const bool fixed_status=false) override
Definition: StandardWell_impl.hpp:2365
StandardWell(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters ¶m, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Definition: StandardWell_impl.hpp:52
void updateWaterThroughput(const double dt, WellState< Scalar > &well_state) const override
Definition: StandardWell_impl.hpp:2025
std::vector< Scalar > computeCurrentWellRates(const Simulator &ebosSimulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2481
void updatePrimaryVariablesNewton(const BVectorWell &dwells, const bool stop_or_zero_rate_target, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:739
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:83
void calculateExplicitQuantities(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1382
void solveEqAndUpdateWellState(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1360
void recoverWellSolutionAndUpdateWellState(const Simulator &simulator, const BVector &x, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1428
void updateWellStateFromPrimaryVariables(WellState< Scalar > &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:762
typename StdWellEval::StdWellConnections::Properties WellConnectionProps
Definition: StandardWell.hpp:262
void computeWellRatesWithBhp(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1448
WellConnectionProps computePropertiesForWellConnectionPressures(const Simulator &simulator, const WellState< Scalar > &well_state) const
Definition: StandardWell_impl.hpp:1121
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:661
void computeWellRatesWithBhpIterations(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1494
void updateIPR(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:781
std::vector< Scalar > computeWellPotentialWithTHP(const Simulator &ebosSimulator, DeferredLogger &deferred_logger, const WellState< Scalar > &well_state) const
Definition: StandardWell_impl.hpp:1558
std::optional< Scalar > computeBhpAtThpLimitProdWithAlq(const Simulator &ebos_simulator, const SummaryState &summary_state, const Scalar alq_value, DeferredLogger &deferred_logger, bool iterate_if_no_solution) const override
Definition: StandardWell_impl.hpp:2204
void computePerfRate(const IntensiveQuantities &intQuants, const std::vector< Value > &mob, const Value &bhp, const std::vector< Scalar > &Tw, const int perf, const bool allow_cf, std::vector< Value > &cq_s, PerforationRates< Scalar > &perf_rates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:94
void updateWellState(const Simulator &simulator, const BVectorWell &dwells, WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:717
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:342
void computeWellConnectionDensitesPressures(const Simulator &simulator, const WellState< Scalar > &well_state, const WellConnectionProps &props, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1291
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1396
std::optional< Scalar > computeBhpAtThpLimitProd(const WellState< Scalar > &well_state, const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2189
void checkConvergenceExtraEqs(const std::vector< Scalar > &res, ConvergenceReport &report) const
Definition: StandardWell_impl.hpp:2127
typename StdWellEval::Eval Eval
Definition: StandardWell.hpp:118
Scalar computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1672
void updateConnectionRatePolyMW(const EvalWell &cq_s_poly, const IntensiveQuantities &int_quants, const WellState< Scalar > &well_state, const int perf, std::vector< RateVector > &connectionRates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2146
bool openCrossFlowAvoidSingularity(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1110
void updatePrimaryVariables(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1793
virtual ConvergenceReport getWellConvergence(const Simulator &simulator, const WellState< Scalar > &well_state, const std::vector< Scalar > &B_avg, DeferredLogger &deferred_logger, const bool relax_tolerance) const override
check whether the well equations get converged for this well
Definition: StandardWell_impl.hpp:1172
std::optional< Scalar > computeBhpAtThpLimitInj(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2286
bool allDrawDownWrongDirection(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1068
EvalWell pskin(const Scalar throughput, const EvalWell &water_velocity, const EvalWell &poly_inj_conc, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1951
static constexpr int numWellConservationEq
Definition: StandardWell.hpp:95
void handleInjectivityEquations(const Simulator &simulator, const WellState< Scalar > &well_state, const int perf, const EvalWell &water_flux_s, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:2079
int setPrimaryVars(typename std::vector< Scalar >::const_iterator it) override
Definition: StandardWell_impl.hpp:2534
void assembleWellEqWithoutIterationImpl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:368
void checkOperabilityUnderTHPLimit(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1019
void addWellPressureEquations(PressureMatrix &mat, const BVector &x, const int pressureVarIndex, const bool use_well_weights, const WellState< Scalar > &well_state) const override
Definition: StandardWell_impl.hpp:1899
EvalWell pskinwater(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1919
void handleInjectivityRate(const Simulator &simulator, const int perf, std::vector< EvalWell > &cq_s) const
Definition: StandardWell_impl.hpp:2056
Scalar getRefDensity() const override
Definition: StandardWell_impl.hpp:1816
void calculateSinglePerf(const Simulator &simulator, const int perf, WellState< Scalar > &well_state, std::vector< RateVector > &connectionRates, std::vector< EvalWell > &cq_s, EvalWell &water_flux_s, EvalWell &cq_s_zfrac_effective, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:494
bool computeWellPotentialsImplicit(const Simulator &ebos_simulator, const WellState< Scalar > &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1593
Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const override
Definition: StandardWell_impl.hpp:1778
EvalWell getQs(const int compIdx) const
Returns scaled rate for a component.
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
std::optional< Scalar > computeBhpAtThpLimitProd(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar maxPerfPress, const Scalar rho, const Scalar alq_value, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for a producer.
std::optional< Scalar > computeBhpAtThpLimitInj(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar rho, const Scalar flo_rel_tol, const int max_iteration, const bool throwOnError, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for an injector.
Scalar calculateThpFromBhp(const std::vector< Scalar > &rates, const Scalar bhp, const Scalar rho, const std::optional< Scalar > &alq, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Calculates THP from BHP.
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
Definition: WellConvergence.hpp:38
void prepareForPotentialCalculations(const SummaryState &summary_state, WellState< Scalar > &well_state, Well::InjectionControls &inj_controls, Well::ProductionControls &prod_controls) const
std::pair< bool, bool > computeWellPotentials(std::vector< Scalar > &well_potentials, const WellState< Scalar > &well_state)
const int num_components_
Definition: WellInterfaceGeneric.hpp:307
void resetDampening()
Definition: WellInterfaceGeneric.hpp:240
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:297
Definition: WellInterfaceIndices.hpp:34
Definition: WellInterface.hpp:77
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:82
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:100
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, 1, 1 > > PressureMatrix
Definition: WellInterface.hpp:97
void getMobility(const Simulator &simulator, const int local_perf_index, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1895
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:86
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:83
bool solveWellWithOperabilityCheck(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:566
Dune::BlockVector< VectorBlockType > BVector
Definition: WellInterface.hpp:96
typename Base::ModelParameters ModelParameters
Definition: WellInterface.hpp:106
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:84
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:85
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: WellInterface.hpp:88
Definition: WellProdIndexCalculator.hpp:37
Scalar connectionProdIndStandard(const std::size_t connIdx, const Scalar connMobility) const
Definition: WellState.hpp:65
const SingleWellState< Scalar > & well(std::size_t well_index) const
Definition: WellState.hpp:285
const PhaseUsage & phaseUsage() const
Definition: WellState.hpp:244
std::vector< Scalar > & wellRates(std::size_t well_index)
One rate per well and phase.
Definition: WellState.hpp:250
@ NONE
Definition: DeferredLogger.hpp:46
Definition: blackoilboundaryratevector.hh:39
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Static data associated with a well perforation.
Definition: PerforationData.hpp:30
Definition: PerforationData.hpp:41
Scalar dis_gas
Definition: PerforationData.hpp:42
Scalar vap_wat
Definition: PerforationData.hpp:45
Scalar vap_oil
Definition: PerforationData.hpp:44
Scalar dis_gas_in_water
Definition: PerforationData.hpp:43
Definition: BlackoilPhases.hpp:46