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_conservation_quantities,
60 const int index_of_well,
62 :
Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_conservation_quantities, num_phases, index_of_well, perf_data)
73 template<
typename TypeTag>
76 init(
const std::vector<Scalar>& depth_arg,
78 const std::vector< Scalar >& B_avg,
79 const bool changed_to_open_this_step)
81 Base::init(depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
82 this->StdWellEval::init(this->perf_depth_, depth_arg, Base::has_polymermw);
89 template<
typename TypeTag>
94 const std::vector<Value>& mob,
96 const std::vector<Scalar>& Tw,
99 std::vector<Value>& cq_s,
103 auto obtain = [
this](
const Eval& value)
105 if constexpr (std::is_same_v<Value, Scalar>) {
106 static_cast<void>(
this);
107 return getValue(value);
109 return this->extendEval(value);
112 auto obtainN = [](
const auto& value)
114 if constexpr (std::is_same_v<Value, Scalar>) {
115 return getValue(value);
120 auto zeroElem = [
this]()
122 if constexpr (std::is_same_v<Value, Scalar>) {
123 static_cast<void>(
this);
126 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
130 const auto& fs = intQuants.fluidState();
131 const Value pressure = obtain(this->getPerfCellPressure(fs));
132 const Value rs = obtain(fs.Rs());
133 const Value rv = obtain(fs.Rv());
134 const Value rvw = obtain(fs.Rvw());
135 const Value rsw = obtain(fs.Rsw());
137 std::vector<Value> b_perfcells_dense(this->numConservationQuantities(), zeroElem());
138 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
139 if (!FluidSystem::phaseIsActive(phaseIdx)) {
142 const unsigned compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
143 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
145 if constexpr (has_solvent) {
146 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
149 if constexpr (has_zFraction) {
150 if (this->isInjector()) {
151 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
152 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
153 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
157 Value skin_pressure = zeroElem();
159 if (this->isInjector()) {
160 const int pskin_index = Bhp + 1 + this->numLocalPerfs() + perf;
161 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
166 std::vector<Value> cmix_s(this->numConservationQuantities(), zeroElem());
167 for (
int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
168 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
191 template<
typename TypeTag>
192 template<
class Value>
196 const Value& pressure,
202 std::vector<Value>& b_perfcells_dense,
203 const std::vector<Scalar>& Tw,
206 const Value& skin_pressure,
207 const std::vector<Value>& cmix_s,
208 std::vector<Value>& cq_s,
213 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
214 Value drawdown = pressure - well_pressure;
215 if (this->isInjector()) {
216 drawdown += skin_pressure;
220 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
221 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx)
223 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
224 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx)
226 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
227 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx)
235 if (!allow_cf && this->isInjector()) {
240 for (
int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
241 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
242 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
245 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
246 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
248 ratioCalc.gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw,
249 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
251 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
252 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
254 ratioCalc.gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw, this->isProducer());
258 if (!allow_cf && this->isProducer()) {
263 Value total_mob_dense = mob[0];
264 for (
int componentIdx = 1; componentIdx < this->numConservationQuantities(); ++componentIdx) {
265 total_mob_dense += mob[componentIdx];
269 Value volumeRatio = bhp * 0.0;
271 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
272 ratioCalc.disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
273 cmix_s, b_perfcells_dense, deferred_logger);
277 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
278 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
279 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
282 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
283 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
284 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
287 if constexpr (Indices::enableSolvent) {
288 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
291 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
292 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
294 ratioCalc.gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
295 cmix_s, b_perfcells_dense,
298 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
299 const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
300 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
302 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
303 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
304 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
310 for (
int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
311 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
312 Value cqt_is = cqt_i / volumeRatio;
313 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
317 if (this->isProducer()) {
318 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
319 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
321 ratioCalc.gasOilPerfRateInj(cq_s, perf_rates,
322 rv, rs, pressure, rvw,
323 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
326 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
327 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
330 ratioCalc.gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
331 pressure, deferred_logger);
338 template<
typename TypeTag>
343 const Well::InjectionControls& inj_controls,
344 const Well::ProductionControls& prod_controls,
351 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
354 this->linSys_.clear();
356 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
357 prod_controls, well_state,
358 group_state, deferred_logger);
364 template<
typename TypeTag>
369 const Well::InjectionControls& inj_controls,
370 const Well::ProductionControls& prod_controls,
376 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
377 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
379 auto& ws = well_state.
well(this->index_of_well_);
380 ws.phase_mixing_rates.fill(0.0);
383 const int np = this->number_of_phases_;
385 std::vector<RateVector> connectionRates = this->connectionRates_;
387 auto& perf_data = ws.perf_data;
388 auto& perf_rates = perf_data.phase_rates;
389 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
391 std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
392 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
393 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
394 calculateSinglePerf(simulator, perf, well_state, connectionRates,
395 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
398 if constexpr (has_polymer && Base::has_polymermw) {
399 if (this->isInjector()) {
400 handleInjectivityEquations(simulator, well_state, perf,
401 water_flux_s, deferred_logger);
404 for (
int componentIdx = 0; componentIdx < this->num_conservation_quantities_; ++componentIdx) {
406 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
408 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
411 assemblePerforationEq(cq_s_effective,
414 this->primary_variables_.numWellEq(),
418 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
419 auto& perf_rate_solvent = perf_data.solvent_rates;
420 perf_rate_solvent[perf] = cq_s[componentIdx].value();
422 perf_rates[perf*np + FluidSystem::activeCompToActivePhaseIdx(componentIdx)] = cq_s[componentIdx].value();
426 if constexpr (has_zFraction) {
428 assembleZFracEq(cq_s_zfrac_effective,
430 this->primary_variables_.numWellEq(),
435 this->connectionRates_ = connectionRates;
440 const auto& comm = this->parallel_well_info_.communication();
441 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
445 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
448 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
451 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
452 if (FluidSystem::numActivePhases() > 1) {
454 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
455 this->F0_[componentIdx]) * volume / dt;
457 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
459 assembleSourceEq(resWell_loc,
461 this->primary_variables_.numWellEq(),
465 const auto& summaryState = simulator.vanguard().summaryState();
466 const Schedule& schedule = simulator.vanguard().schedule();
467 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
469 assembleControlEq(well_state, group_state,
470 schedule, summaryState,
471 inj_controls, prod_controls,
472 this->primary_variables_,
473 this->connections_.rho(),
475 stopped_or_zero_target,
481 this->linSys_.invert();
483 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
490 template<
typename TypeTag>
496 std::vector<RateVector>& connectionRates,
497 std::vector<EvalWell>& cq_s,
502 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
503 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
504 const int cell_idx = this->well_cells_[perf];
505 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
506 std::vector<EvalWell> mob(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
507 getMobility(simulator, perf, mob, deferred_logger);
510 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
511 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
512 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
513 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
514 cq_s, perf_rates, deferred_logger);
516 auto& ws = well_state.
well(this->index_of_well_);
517 auto& perf_data = ws.perf_data;
518 if constexpr (has_polymer && Base::has_polymermw) {
519 if (this->isInjector()) {
522 const unsigned water_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
523 water_flux_s = cq_s[water_comp_idx];
526 handleInjectivityRate(simulator, perf, cq_s);
531 if (this->isProducer()) {
532 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.
dis_gas;
533 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.
dis_gas_in_water;
534 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.
vap_oil;
535 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.
vap_wat;
536 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.
dis_gas;
537 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.
dis_gas_in_water;
538 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.
vap_oil;
539 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.
vap_wat;
542 if constexpr (has_energy) {
543 connectionRates[perf][Indices::contiEnergyEqIdx] =
544 connectionRateEnergy(cq_s, intQuants, deferred_logger);
547 if constexpr (has_polymer) {
548 std::variant<Scalar,EvalWell> polymerConcentration;
549 if (this->isInjector()) {
550 polymerConcentration = this->wpolymer();
552 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
553 intQuants.polymerViscosityCorrection());
556 [[maybe_unused]]
EvalWell cq_s_poly;
557 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
559 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
560 cq_s, polymerConcentration);
562 if constexpr (Base::has_polymermw) {
563 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
564 perf, connectionRates, deferred_logger);
568 if constexpr (has_foam) {
569 std::variant<Scalar,EvalWell> foamConcentration;
570 if (this->isInjector()) {
571 foamConcentration = this->wfoam();
573 foamConcentration = this->extendEval(intQuants.foamConcentration());
575 connectionRates[perf][Indices::contiFoamEqIdx] =
576 this->connections_.connectionRateFoam(cq_s, foamConcentration,
577 FoamModule::transportPhase(),
581 if constexpr (has_zFraction) {
582 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
583 if (this->isInjector()) {
584 solventConcentration = this->wsolvent();
586 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
587 this->extendEval(intQuants.yVolume())};
589 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
590 cq_s_zfrac_effective) =
591 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
593 solventConcentration);
596 if constexpr (has_brine) {
597 std::variant<Scalar,EvalWell> saltConcentration;
598 if (this->isInjector()) {
599 saltConcentration = this->wsalt();
601 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
604 connectionRates[perf][Indices::contiBrineEqIdx] =
605 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
610 if constexpr (has_micp) {
611 std::variant<Scalar,EvalWell> microbialConcentration;
612 std::variant<Scalar,EvalWell> oxygenConcentration;
613 std::variant<Scalar,EvalWell> ureaConcentration;
614 if (this->isInjector()) {
615 microbialConcentration = this->wmicrobes();
616 oxygenConcentration = this->woxygen();
617 ureaConcentration = this->wurea();
619 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
620 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
621 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
623 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
624 connectionRates[perf][Indices::contiOxygenEqIdx],
625 connectionRates[perf][Indices::contiUreaEqIdx]) =
626 this->connections_.connectionRatesMICP(perf_data.microbial_rates[perf],
627 perf_data.oxygen_rates[perf],
628 perf_data.urea_rates[perf],
630 microbialConcentration,
636 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
639 if (FluidSystem::phaseUsage().hasCO2orH2Store()) {
640 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
641 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
642 perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
646 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
647 const unsigned wat_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
648 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::waterPhaseIdx, Base::pvtRegionIdx() );
649 perf_data.wat_mass_rates[perf] = cq_s[wat_comp_idx].value() * rho;
655 template<
typename TypeTag>
656 template<
class Value>
661 std::vector<Value>& mob,
664 auto obtain = [
this](
const Eval& value)
666 if constexpr (std::is_same_v<Value, Scalar>) {
667 static_cast<void>(
this);
668 return getValue(value);
670 return this->extendEval(value);
674 obtain, deferred_logger);
677 if constexpr (has_polymer) {
678 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
679 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
684 if constexpr (!Base::has_polymermw) {
685 if constexpr (std::is_same_v<Value, Scalar>) {
686 std::vector<EvalWell> mob_eval(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
687 for (std::size_t i = 0; i < mob.size(); ++i) {
688 mob_eval[i].setValue(mob[i]);
690 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
691 for (std::size_t i = 0; i < mob.size(); ++i) {
692 mob[i] = getValue(mob_eval[i]);
695 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
702 const Scalar bhp = this->primary_variables_.value(Bhp);
703 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
704 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
705 for (std::size_t i = 0; i < mob.size(); ++i) {
706 mob[i] *= multiplier;
712 template<
typename TypeTag>
720 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
722 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
723 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
725 const auto& summary_state = simulator.vanguard().summaryState();
726 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
727 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.
well(this->index_of_well_));
734 template<
typename TypeTag>
738 const bool stop_or_zero_rate_target,
741 const Scalar dFLimit = this->param_.dwell_fraction_max_;
742 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
743 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
746 if constexpr (Base::has_polymermw) {
747 this->primary_variables_.updateNewtonPolyMW(dwells);
750 this->primary_variables_.checkFinite(deferred_logger);
757 template<
typename TypeTag>
761 const SummaryState& summary_state,
764 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
767 if constexpr (Base::has_polymermw) {
768 this->primary_variables_.copyToWellStatePolyMW(well_state);
776 template<
typename TypeTag>
784 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
785 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
787 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
788 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
789 getMobility(simulator, perf, mob, deferred_logger);
791 const int cell_idx = this->well_cells_[perf];
792 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
793 const auto& fs = int_quantities.fluidState();
795 Scalar p_r = this->getPerfCellPressure(fs).value();
798 std::vector<Scalar> b_perf(this->num_conservation_quantities_);
799 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
800 if (!FluidSystem::phaseIsActive(phase)) {
803 const unsigned comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phase));
804 b_perf[comp_idx] = fs.invB(phase).value();
806 if constexpr (has_solvent) {
807 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
811 const Scalar h_perf = this->connections_.pressure_diff(perf);
812 const Scalar pressure_diff = p_r - h_perf;
817 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
818 deferred_logger.
debug(
"CROSSFLOW_IPR",
819 "cross flow found when updateIPR for well " + name()
820 +
" . The connection is ignored in IPR calculations");
826 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
827 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
828 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
832 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
833 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
834 for (
int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
835 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
836 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
837 ipr_b_perf[comp_idx] += tw_mob;
841 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
842 const unsigned oil_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
843 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
844 const Scalar rs = (fs.Rs()).value();
845 const Scalar rv = (fs.Rv()).value();
847 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
848 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
850 ipr_a_perf[gas_comp_idx] += dis_gas_a;
851 ipr_a_perf[oil_comp_idx] += vap_oil_a;
853 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
854 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
856 ipr_b_perf[gas_comp_idx] += dis_gas_b;
857 ipr_b_perf[oil_comp_idx] += vap_oil_b;
860 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
861 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
862 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
865 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
866 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
869 template<
typename TypeTag>
883 auto rates = well_state.
well(this->index_of_well_).surface_rates;
885 for (std::size_t p = 0; p < rates.size(); ++p) {
886 zero_rates &= rates[p] == 0.0;
888 auto& ws = well_state.
well(this->index_of_well_);
890 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
891 deferred_logger.
debug(msg);
903 const auto& group_state = simulator.problem().wellModel().groupState();
905 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
906 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
908 auto inj_controls = Well::InjectionControls(0);
909 auto prod_controls = Well::ProductionControls(0);
910 prod_controls.addControl(Well::ProducerCMode::BHP);
911 prod_controls.bhp_limit = well_state.
well(this->index_of_well_).bhp;
914 const auto cmode = ws.production_cmode;
915 ws.production_cmode = Well::ProducerCMode::BHP;
916 const double dt = simulator.timeStepSize();
917 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
919 const size_t nEq = this->primary_variables_.numWellEq();
923 for (
size_t i=0; i < nEq; ++i){
929 x_well[0].resize(nEq);
930 this->linSys_.solve(rhs, x_well);
932 for (
int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
933 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
934 const int idx = FluidSystem::activeCompToActivePhaseIdx(comp_idx);
935 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
937 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
939 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
942 ws.production_cmode = cmode;
945 template<
typename TypeTag>
952 const auto& summaryState = simulator.vanguard().summaryState();
956 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
957 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
960 Scalar total_ipr_mass_rate = 0.0;
961 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
963 if (!FluidSystem::phaseIsActive(phaseIdx)) {
967 const unsigned compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
968 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
970 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
971 total_ipr_mass_rate += ipr_rate * rho;
973 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
974 this->operability_status_.operable_under_only_bhp_limit =
false;
978 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
982 std::vector<Scalar> well_rates_bhp_limit;
983 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
985 this->adaptRatesForVFP(well_rates_bhp_limit);
986 const Scalar thp_limit = this->getTHPConstraint(summaryState);
989 this->connections_.rho(),
990 this->getALQ(well_state),
993 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
994 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1005 this->operability_status_.operable_under_only_bhp_limit =
true;
1006 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1014 template<
typename TypeTag>
1021 const auto& summaryState = simulator.vanguard().summaryState();
1022 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1023 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1026 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1029 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1030 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1032 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1033 if (this->isProducer() && *obtain_bhp < thp_limit) {
1034 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1035 +
" bars is SMALLER than thp limit "
1037 +
" bars as a producer for well " + name();
1038 deferred_logger.
debug(msg);
1040 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1041 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1042 +
" bars is LARGER than thp limit "
1044 +
" bars as a injector for well " + name();
1045 deferred_logger.
debug(msg);
1048 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1049 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1050 if (!this->wellIsStopped()) {
1051 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1052 deferred_logger.
debug(
" could not find bhp value at thp limit "
1054 +
" bar for well " + name() +
", the well might need to be closed ");
1063 template<
typename TypeTag>
1068 bool all_drawdown_wrong_direction =
true;
1070 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1071 const int cell_idx = this->well_cells_[perf];
1072 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1073 const auto& fs = intQuants.fluidState();
1075 const Scalar pressure = this->getPerfCellPressure(fs).value();
1076 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1079 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1080 const Scalar drawdown = pressure - well_pressure;
1085 if ( (drawdown < 0. && this->isInjector()) ||
1086 (drawdown > 0. && this->isProducer()) ) {
1087 all_drawdown_wrong_direction =
false;
1092 const auto& comm = this->parallel_well_info_.communication();
1093 if (comm.size() > 1)
1095 all_drawdown_wrong_direction =
1096 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1099 return all_drawdown_wrong_direction;
1105 template<
typename TypeTag>
1110 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1116 template<
typename TypeTag>
1122 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1124 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1126 return model.intensiveQuantities(cell_idx, 0)
1127 .fluidState().temperature(phase_idx).value();
1131 [&model = simulator.model()](
int cell_idx)
1133 return model.intensiveQuantities(cell_idx, 0)
1134 .fluidState().saltConcentration().value();
1138 [&model = simulator.model()](
int cell_idx)
1140 return model.intensiveQuantities(cell_idx, 0)
1141 .fluidState().pvtRegionIndex();
1145 if constexpr (Indices::enableSolvent) {
1146 prop_func.solventInverseFormationVolumeFactor =
1147 [&model = simulator.model()](
int cell_idx)
1149 return model.intensiveQuantities(cell_idx, 0)
1150 .solventInverseFormationVolumeFactor().value();
1153 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1155 return model.intensiveQuantities(cell_idx, 0)
1156 .solventRefDensity();
1160 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1167 template<
typename TypeTag>
1172 const std::vector<Scalar>& B_avg,
1174 const bool relax_tolerance)
const
1178 assert((
int(B_avg.size()) == this->num_conservation_quantities_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1180 Scalar tol_wells = this->param_.tolerance_wells_;
1182 constexpr Scalar stopped_factor = 1.e-4;
1184 constexpr Scalar dynamic_thp_factor = 1.e-1;
1185 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1186 tol_wells = tol_wells*stopped_factor;
1187 }
else if (this->getDynamicThpLimit()) {
1188 tol_wells = tol_wells*dynamic_thp_factor;
1191 std::vector<Scalar> res;
1194 this->param_.max_residual_allowed_,
1196 this->param_.relaxed_tolerance_flow_well_,
1198 this->wellIsStopped(),
1202 checkConvergenceExtraEqs(res, report);
1211 template<
typename TypeTag>
1219 auto fluidState = [&simulator,
this](
const int perf)
1221 const auto cell_idx = this->well_cells_[perf];
1222 return simulator.model()
1223 .intensiveQuantities(cell_idx, 0).fluidState();
1226 const int np = this->number_of_phases_;
1227 auto setToZero = [np](
Scalar* x) ->
void
1229 std::fill_n(x, np, 0.0);
1232 auto addVector = [np](
const Scalar* src,
Scalar* dest) ->
void
1234 std::transform(src, src + np, dest, dest, std::plus<>{});
1237 auto& ws = well_state.
well(this->index_of_well_);
1238 auto& perf_data = ws.perf_data;
1239 auto* wellPI = ws.productivity_index.data();
1240 auto* connPI = perf_data.prod_index.data();
1244 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1245 auto subsetPerfID = 0;
1247 for (
const auto& perf : *this->perf_data_) {
1248 auto allPerfID = perf.ecl_index;
1250 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) ->
Scalar
1255 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
1256 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1258 const auto& fs = fluidState(subsetPerfID);
1261 if (this->isInjector()) {
1262 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1263 mob, connPI, deferred_logger);
1266 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1269 addVector(connPI, wellPI);
1276 const auto& comm = this->parallel_well_info_.communication();
1277 if (comm.size() > 1) {
1278 comm.sum(wellPI, np);
1281 assert ((
static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1282 "Internal logic error in processing connections for PI/II");
1287 template<
typename TypeTag>
1297 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1302 [&model = simulator.model()](
const int cell,
1303 const std::vector<int>& phases,
1304 std::vector<Scalar>& mob)
1306 const auto& iq = model.intensiveQuantities(cell, 0);
1308 std::transform(phases.begin(), phases.end(), mob.begin(),
1309 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1314 [&model = simulator.model()](
const int cell,
1315 const std::vector<int>& phases,
1316 std::vector<Scalar>& rho)
1318 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1320 std::transform(phases.begin(), phases.end(), rho.begin(),
1321 [&fs](
const int phase) { return fs.density(phase).value(); });
1325 const auto stopped_or_zero_rate_target = this->
1326 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1329 .computeProperties(stopped_or_zero_rate_target, well_state,
1330 prop_func, props, deferred_logger);
1337 template<
typename TypeTag>
1344 const auto props = computePropertiesForWellConnectionPressures
1345 (simulator, well_state);
1347 computeWellConnectionDensitesPressures(simulator, well_state,
1348 props, deferred_logger);
1355 template<
typename TypeTag>
1362 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1367 dx_well[0].resize(this->primary_variables_.numWellEq());
1368 this->linSys_.solve( dx_well);
1370 updateWellState(simulator, dx_well, well_state, deferred_logger);
1377 template<
typename TypeTag>
1384 updatePrimaryVariables(simulator, well_state, deferred_logger);
1385 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1386 this->computeAccumWell();
1391 template<
typename TypeTag>
1396 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1398 if (this->param_.matrix_add_well_contributions_)
1404 this->linSys_.apply(x, Ax);
1410 template<
typename TypeTag>
1415 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1417 this->linSys_.apply(r);
1423 template<
typename TypeTag>
1431 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1434 xw[0].resize(this->primary_variables_.numWellEq());
1436 this->linSys_.recoverSolutionWell(x, xw);
1437 updateWellState(simulator, xw, well_state, deferred_logger);
1443 template<
typename TypeTag>
1448 std::vector<Scalar>& well_flux,
1452 const int np = this->number_of_phases_;
1453 well_flux.resize(np, 0.0);
1455 const bool allow_cf = this->getAllowCrossFlow();
1457 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1458 const int cell_idx = this->well_cells_[perf];
1459 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1461 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
1462 getMobility(simulator, perf, mob, deferred_logger);
1463 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1464 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1465 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1467 std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
1469 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1470 cq_s, perf_rates, deferred_logger);
1472 for(
int p = 0; p < np; ++p) {
1473 well_flux[FluidSystem::activeCompToActivePhaseIdx(p)] += cq_s[p];
1477 if constexpr (has_solvent) {
1478 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
1480 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1481 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1484 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1489 template<
typename TypeTag>
1494 std::vector<Scalar>& well_flux,
1505 WellStateType well_state_copy = simulator.problem().wellModel().wellState();
1506 const auto& group_state = simulator.problem().wellModel().groupState();
1509 const auto& summary_state = simulator.vanguard().summaryState();
1510 auto inj_controls = well_copy.
well_ecl_.isInjector()
1511 ? well_copy.
well_ecl_.injectionControls(summary_state)
1512 : Well::InjectionControls(0);
1513 auto prod_controls = well_copy.
well_ecl_.isProducer()
1514 ? well_copy.
well_ecl_.productionControls(summary_state) :
1515 Well::ProductionControls(0);
1518 auto& ws = well_state_copy.
well(this->index_of_well_);
1520 inj_controls.bhp_limit = bhp;
1521 ws.injection_cmode = Well::InjectorCMode::BHP;
1523 prod_controls.bhp_limit = bhp;
1524 ws.production_cmode = Well::ProducerCMode::BHP;
1529 const int np = this->number_of_phases_;
1530 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1531 for (
int phase = 0; phase < np; ++phase){
1532 well_state_copy.
wellRates(this->index_of_well_)[phase]
1533 = sign * ws.well_potentials[phase];
1538 const double dt = simulator.timeStepSize();
1539 const bool converged = well_copy.
iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1541 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1542 " potentials are computed based on unconverged solution";
1543 deferred_logger.
debug(msg);
1553 template<
typename TypeTag>
1554 std::vector<typename StandardWell<TypeTag>::Scalar>
1560 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1561 const auto& summary_state = simulator.vanguard().summaryState();
1563 const auto& well = this->well_ecl_;
1564 if (well.isInjector()){
1565 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1566 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1567 if (bhp_at_thp_limit) {
1568 const Scalar bhp = std::min(*bhp_at_thp_limit,
1569 static_cast<Scalar>(controls.bhp_limit));
1570 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1572 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1573 "Failed in getting converged thp based potential calculation for well "
1574 + name() +
". Instead the bhp based value is used");
1575 const Scalar bhp = controls.bhp_limit;
1576 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1579 computeWellRatesWithThpAlqProd(
1580 simulator, summary_state,
1581 deferred_logger, potentials, this->getALQ(well_state)
1588 template<
typename TypeTag>
1593 std::vector<Scalar>& well_potentials,
1603 const auto& group_state = simulator.problem().wellModel().groupState();
1604 auto& ws = well_state_copy.
well(this->index_of_well_);
1607 const auto& summary_state = simulator.vanguard().summaryState();
1608 auto inj_controls = well_copy.
well_ecl_.isInjector()
1609 ? well_copy.
well_ecl_.injectionControls(summary_state)
1610 : Well::InjectionControls(0);
1611 auto prod_controls = well_copy.
well_ecl_.isProducer()
1612 ? well_copy.
well_ecl_.productionControls(summary_state) :
1613 Well::ProductionControls(0);
1619 const int num_perf = ws.perf_data.size();
1620 for (
int perf = 0; perf < num_perf; ++perf) {
1624 const int np = this->number_of_phases_;
1625 bool trivial =
true;
1626 for (
int phase = 0; phase < np; ++phase){
1627 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1631 for (
int phase = 0; phase < np; ++phase) {
1632 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1637 const double dt = simulator.timeStepSize();
1639 bool converged =
false;
1640 if (this->well_ecl_.isProducer()) {
1641 converged = well_copy.
solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1643 converged = well_copy.
iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1647 well_potentials.clear();
1648 well_potentials.resize(np, 0.0);
1649 for (
int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
1650 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1652 well_potentials[FluidSystem::activeCompToActivePhaseIdx(comp_idx)] = rate.value();
1656 if constexpr (has_solvent) {
1657 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
1659 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1661 well_potentials[gas_pos] += rate.value();
1667 template<
typename TypeTag>
1671 const SummaryState &summary_state,
1673 std::vector<Scalar>& potentials,
1677 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1678 simulator, summary_state, alq, deferred_logger,
true);
1679 if (bhp_at_thp_limit) {
1680 const auto& controls = this->well_ecl_.productionControls(summary_state);
1681 bhp = std::max(*bhp_at_thp_limit,
1682 static_cast<Scalar>(controls.bhp_limit));
1683 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1686 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1687 "Failed in getting converged thp based potential calculation for well "
1688 + name() +
". Instead the bhp based value is used");
1689 const auto& controls = this->well_ecl_.productionControls(summary_state);
1690 bhp = controls.bhp_limit;
1691 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1696 template<
typename TypeTag>
1700 const SummaryState& summary_state,
1702 std::vector<Scalar>& potentials,
1706 computeWellRatesAndBhpWithThpAlqProd(simulator,
1713 template<
typename TypeTag>
1718 std::vector<Scalar>& well_potentials,
1721 const auto [compute_potential, bhp_controlled_well] =
1724 if (!compute_potential) {
1728 bool converged_implicit =
false;
1732 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1733 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1735 if (!converged_implicit) {
1737 const auto& summaryState = simulator.vanguard().summaryState();
1738 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1748 const auto& ws = well_state.
well(this->index_of_well_);
1749 if (this->isInjector())
1750 bhp = std::max(ws.bhp, bhp);
1752 bhp = std::min(ws.bhp, bhp);
1754 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1755 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1758 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1762 this->checkNegativeWellPotentials(well_potentials,
1763 this->param_.check_well_operability_,
1773 template<
typename TypeTag>
1777 const int openConnIdx)
const
1779 return (openConnIdx < 0)
1781 : this->connections_.rho(openConnIdx);
1788 template<
typename TypeTag>
1795 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1797 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1798 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1801 if constexpr (Base::has_polymermw) {
1802 this->primary_variables_.updatePolyMW(well_state);
1805 this->primary_variables_.checkFinite(deferred_logger);
1811 template<
typename TypeTag>
1816 return this->connections_.rho();
1822 template<
typename TypeTag>
1827 std::vector<EvalWell>& mob,
1830 const int cell_idx = this->well_cells_[perf];
1831 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1832 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1836 if (this->isInjector()) {
1838 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1839 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
1840 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1843 if (PolymerModule::hasPlyshlog()) {
1846 if (this->isInjector() && this->wpolymer() == 0.) {
1851 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1852 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1854 std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1856 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1857 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1858 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1859 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1860 perf_rates, deferred_logger);
1862 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1863 const auto& material_law_manager = simulator.problem().materialLawManager();
1864 const auto& scaled_drainage_info =
1865 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1866 const Scalar swcr = scaled_drainage_info.Swcr;
1867 const EvalWell poro = this->extendEval(int_quant.porosity());
1868 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1870 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1871 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
1872 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1874 if (PolymerModule::hasShrate()) {
1877 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1879 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1880 int_quant.pvtRegionIndex(),
1883 mob[waterCompIdx] /= shear_factor;
1887 template<
typename TypeTag>
1891 this->linSys_.extract(jacobian);
1895 template <
typename TypeTag>
1899 const int pressureVarIndex,
1900 const bool use_well_weights,
1903 this->linSys_.extractCPRPressureMatrix(jacobian,
1914 template<
typename TypeTag>
1921 if constexpr (Base::has_polymermw) {
1922 const int water_table_id = this->polymerWaterTable_();
1923 if (water_table_id <= 0) {
1925 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1928 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1929 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1931 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1932 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1936 fmt::format(
"Polymermw is not activated, while injecting "
1937 "skin pressure is requested for well {}", name()),
1946 template<
typename TypeTag>
1954 if constexpr (Base::has_polymermw) {
1955 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1956 const EvalWell water_velocity_abs = abs(water_velocity);
1957 if (poly_inj_conc == 0.) {
1958 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1960 const int polymer_table_id = this->polymerTable_();
1961 if (polymer_table_id <= 0) {
1963 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1966 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1967 const Scalar reference_concentration = skprpolytable.refConcentration;
1968 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1970 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1971 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1972 if (poly_inj_conc == reference_concentration) {
1973 return sign * pskin_poly;
1976 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1977 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1978 return sign * pskin;
1981 fmt::format(
"Polymermw is not activated, while injecting "
1982 "skin pressure is requested for well {}", name()),
1991 template<
typename TypeTag>
1998 if constexpr (Base::has_polymermw) {
1999 const int table_id = this->polymerInjTable_();
2000 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2001 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2002 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2003 if (this->wpolymer() == 0.) {
2004 return molecular_weight;
2006 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2007 return molecular_weight;
2010 fmt::format(
"Polymermw is not activated, while injecting "
2011 "polymer molecular weight is requested for well {}", name()),
2020 template<
typename TypeTag>
2026 if constexpr (Base::has_polymermw) {
2027 if (!this->isInjector()) {
2031 auto& perf_water_throughput = well_state.
well(this->index_of_well_)
2032 .perf_data.water_throughput;
2034 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2035 const Scalar perf_water_vel =
2036 this->primary_variables_.value(Bhp + 1 + perf);
2040 if (perf_water_vel >
Scalar{0}) {
2041 perf_water_throughput[perf] += perf_water_vel * dt;
2051 template<
typename TypeTag>
2056 std::vector<EvalWell>& cq_s)
const
2058 const int cell_idx = this->well_cells_[perf];
2059 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2060 const auto& fs = int_quants.fluidState();
2061 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2062 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2063 const int wat_vel_index = Bhp + 1 + perf;
2064 const unsigned water_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
2068 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2074 template<
typename TypeTag>
2083 const int cell_idx = this->well_cells_[perf];
2084 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2085 const auto& fs = int_quants.fluidState();
2086 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2087 const EvalWell water_flux_r = water_flux_s / b_w;
2088 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2089 const EvalWell water_velocity = water_flux_r / area;
2090 const int wat_vel_index = Bhp + 1 + perf;
2093 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2095 const auto& ws = well_state.
well(this->index_of_well_);
2096 const auto& perf_data = ws.perf_data;
2097 const auto& perf_water_throughput = perf_data.water_throughput;
2098 const Scalar throughput = perf_water_throughput[perf];
2099 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2101 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2102 poly_conc.setValue(this->wpolymer());
2105 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2106 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2109 assembleInjectivityEq(eq_pskin,
2114 this->primary_variables_.numWellEq(),
2122 template<
typename TypeTag>
2131 if constexpr (Base::has_polymermw) {
2133 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2141 template<
typename TypeTag>
2148 std::vector<RateVector>& connectionRates,
2153 if (this->isInjector()) {
2154 const int wat_vel_index = Bhp + 1 + perf;
2155 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2156 if (water_velocity > 0.) {
2157 const auto& ws = well_state.
well(this->index_of_well_);
2158 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2159 const Scalar throughput = perf_water_throughput[perf];
2160 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2161 cq_s_polymw *= molecular_weight;
2167 }
else if (this->isProducer()) {
2168 if (cq_s_polymw < 0.) {
2169 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2176 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2183 template<
typename TypeTag>
2184 std::optional<typename StandardWell<TypeTag>::Scalar>
2188 const SummaryState& summary_state,
2191 return computeBhpAtThpLimitProdWithAlq(simulator,
2193 this->getALQ(well_state),
2198 template<
typename TypeTag>
2199 std::optional<typename StandardWell<TypeTag>::Scalar>
2202 const SummaryState& summary_state,
2205 bool iterate_if_no_solution)
const
2209 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2215 std::vector<Scalar> rates(3);
2216 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2217 this->adaptRatesForVFP(rates);
2221 Scalar max_pressure = 0.0;
2222 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2223 const int cell_idx = this->well_cells_[perf];
2224 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2225 const auto& fs = int_quants.fluidState();
2226 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2227 max_pressure = std::max(max_pressure, pressure_cell);
2232 this->connections_.rho(),
2234 this->getTHPConstraint(summary_state),
2238 auto v = frates(*bhpAtLimit);
2239 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2244 if (!iterate_if_no_solution)
2245 return std::nullopt;
2247 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2251 std::vector<Scalar> rates(3);
2252 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2253 this->adaptRatesForVFP(rates);
2260 this->connections_.rho(),
2262 this->getTHPConstraint(summary_state),
2268 auto v = frates(*bhpAtLimit);
2269 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2275 return std::nullopt;
2280 template<
typename TypeTag>
2281 std::optional<typename StandardWell<TypeTag>::Scalar>
2284 const SummaryState& summary_state,
2288 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2294 std::vector<Scalar> rates(3);
2295 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2301 this->connections_.rho(),
2312 template<
typename TypeTag>
2317 const Well::InjectionControls& inj_controls,
2318 const Well::ProductionControls& prod_controls,
2323 updatePrimaryVariables(simulator, well_state, deferred_logger);
2325 const int max_iter = this->param_.max_inner_iter_wells_;
2328 bool relax_convergence =
false;
2329 this->regularize_ =
false;
2331 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2333 if (it > this->param_.strict_inner_iter_wells_) {
2334 relax_convergence =
true;
2335 this->regularize_ =
true;
2338 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2340 converged = report.converged();
2346 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2353 }
while (it < max_iter);
2359 template<
typename TypeTag>
2364 const Well::InjectionControls& inj_controls,
2365 const Well::ProductionControls& prod_controls,
2369 const bool fixed_control ,
2370 const bool fixed_status )
2372 updatePrimaryVariables(simulator, well_state, deferred_logger);
2374 const int max_iter = this->param_.max_inner_iter_wells_;
2376 bool converged =
false;
2377 bool relax_convergence =
false;
2378 this->regularize_ =
false;
2379 const auto& summary_state = simulator.vanguard().summaryState();
2384 constexpr int min_its_after_switch = 4;
2386 const int max_status_switch = this->param_.max_well_status_switch_;
2387 int its_since_last_switch = min_its_after_switch;
2388 int switch_count= 0;
2390 const auto well_status_orig = this->wellStatus_;
2391 const auto operability_orig = this->operability_status_;
2392 auto well_status_cur = well_status_orig;
2393 int status_switch_count = 0;
2395 const bool allow_open = well_state.
well(this->index_of_well_).status == WellStatus::OPEN;
2397 const bool allow_switching =
2398 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2399 (!fixed_control || !fixed_status) && allow_open;
2401 bool changed =
false;
2402 bool final_check =
false;
2404 this->operability_status_.resetOperability();
2405 this->operability_status_.solvable =
true;
2407 its_since_last_switch++;
2408 if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
2409 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2410 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2411 inj_controls, prod_controls, wqTotal,
2412 deferred_logger, fixed_control, fixed_status);
2414 its_since_last_switch = 0;
2416 if (well_status_cur != this->wellStatus_) {
2417 well_status_cur = this->wellStatus_;
2418 status_switch_count++;
2421 if (!changed && final_check) {
2424 final_check =
false;
2426 if (status_switch_count == max_status_switch) {
2427 this->wellStatus_ = well_status_orig;
2431 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2433 if (it > this->param_.strict_inner_iter_wells_) {
2434 relax_convergence =
true;
2435 this->regularize_ =
true;
2438 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2440 converged = report.converged();
2444 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2446 its_since_last_switch = min_its_after_switch;
2453 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2455 }
while (it < max_iter);
2458 if (allow_switching){
2460 const bool is_stopped = this->wellIsStopped();
2461 if (this->wellHasTHPConstraints(summary_state)){
2462 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2463 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2465 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2469 this->wellStatus_ = well_status_orig;
2470 this->operability_status_ = operability_orig;
2471 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2472 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2473 deferred_logger.
debug(message);
2479 template<
typename TypeTag>
2480 std::vector<typename StandardWell<TypeTag>::Scalar>
2486 std::vector<Scalar> well_q_s(this->num_conservation_quantities_, 0.);
2487 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2488 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2489 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2490 const int cell_idx = this->well_cells_[perf];
2491 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2492 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
2493 getMobility(simulator, perf, mob, deferred_logger);
2494 std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
2495 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2496 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2497 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2499 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2500 cq_s, perf_rates, deferred_logger);
2501 for (
int comp = 0; comp < this->num_conservation_quantities_; ++comp) {
2502 well_q_s[comp] += cq_s[comp];
2505 const auto& comm = this->parallel_well_info_.communication();
2506 if (comm.size() > 1)
2508 comm.sum(well_q_s.data(), well_q_s.size());
2515 template <
typename TypeTag>
2516 std::vector<typename StandardWell<TypeTag>::Scalar>
2520 const int num_pri_vars = this->primary_variables_.numWellEq();
2521 std::vector<Scalar> retval(num_pri_vars);
2522 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2523 retval[ii] = this->primary_variables_.value(ii);
2532 template <
typename TypeTag>
2537 const int num_pri_vars = this->primary_variables_.numWellEq();
2538 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2539 this->primary_variables_.setValue(ii, it[ii]);
2541 return num_pri_vars;
2545 template <
typename TypeTag>
2549 const IntensiveQuantities& intQuants,
2552 auto fs = intQuants.fluidState();
2554 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2555 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2560 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2561 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
2562 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2563 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2564 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2567 const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
2568 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
2573 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2575 deferred_logger.
debug(
2576 fmt::format(
"Problematic d value {} obtained for well {}"
2577 " during calculateSinglePerf with rs {}"
2578 ", rv {}. Continue as if no dissolution (rs = 0) and"
2579 " vaporization (rv = 0) for this connection.",
2580 d, this->name(), fs.Rs(), fs.Rv()));
2581 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2583 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2584 cq_r_thermal = (cq_s[gasCompIdx] -
2585 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2586 (d * this->extendEval(fs.invB(phaseIdx)) );
2587 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2589 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2591 (d * this->extendEval(fs.invB(phaseIdx)) );
2597 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2599 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2600 fs.setTemperature(this->well_ecl_.inj_temperature());
2601 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2602 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2603 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2604 paramCache.setRegionIndex(pvtRegionIdx);
2605 paramCache.updatePhase(fs, phaseIdx);
2607 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2608 fs.setDensity(phaseIdx, rho);
2609 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2610 fs.setEnthalpy(phaseIdx, h);
2611 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2612 result += getValue(cq_r_thermal);
2613 }
else if (cq_r_thermal > 0.0) {
2614 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2615 result += Base::restrictEval(cq_r_thermal);
2618 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2619 result += Base::restrictEval(cq_r_thermal);
2623 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:41
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:198
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:106
StdWellConnections connections_
Connection level values.
Definition: StandardWellEval.hpp:105
PrimaryVariables primary_variables_
Primary variables for well.
Definition: StandardWellEval.hpp:99
Definition: StandardWell.hpp:60
void calculateExplicitQuantities(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1380
EvalWell wpolymermw(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1994
typename StdWellEval::EvalWell EvalWell
Definition: StandardWell.hpp:119
void updateWellStateFromPrimaryVariables(WellStateType &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:760
WellConnectionProps computePropertiesForWellConnectionPressures(const Simulator &simulator, const WellStateType &well_state) const
Definition: StandardWell_impl.hpp:1119
typename StdWellEval::BVectorWell BVectorWell
Definition: StandardWell.hpp:120
std::optional< Scalar > computeBhpAtThpLimitProd(const WellStateType &well_state, const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2186
void computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1699
void addWellContributions(SparseMatrixAdapter &mat) const override
Definition: StandardWell_impl.hpp:1889
std::vector< Scalar > getPrimaryVars() const override
Definition: StandardWell_impl.hpp:2518
void addWellPressureEquations(PressureMatrix &mat, const BVector &x, const int pressureVarIndex, const bool use_well_weights, const WellStateType &well_state) const override
Definition: StandardWell_impl.hpp:1897
void updateWaterMobilityWithPolymer(const Simulator &simulator, const int perf, std::vector< EvalWell > &mob_water, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1825
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:341
std::vector< Scalar > computeCurrentWellRates(const Simulator &ebosSimulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2482
void calculateSinglePerf(const Simulator &simulator, const int perf, WellStateType &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:493
void updatePrimaryVariablesNewton(const BVectorWell &dwells, const bool stop_or_zero_rate_target, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:737
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:82
void computeWellConnectionPressures(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1340
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_conservation_quantities, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Definition: StandardWell_impl.hpp:52
typename StdWellEval::StdWellConnections::Properties WellConnectionProps
Definition: StandardWell.hpp:263
void updateConnectionRatePolyMW(const EvalWell &cq_s_poly, const IntensiveQuantities &int_quants, const WellStateType &well_state, const int perf, std::vector< RateVector > &connectionRates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2144
bool iterateWellEqWithSwitching(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &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:2362
void computeWellRatesWithBhp(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1446
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:659
void computeWellRatesWithBhpIterations(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1492
void updateIPR(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:779
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:2201
void updateIPRImplicit(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:872
void computeWellConnectionDensitesPressures(const Simulator &simulator, const WellStateType &well_state, const WellConnectionProps &props, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1289
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:93
void updateWellState(const Simulator &simulator, const BVectorWell &dwells, WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:715
void handleInjectivityEquations(const Simulator &simulator, const WellStateType &well_state, const int perf, const EvalWell &water_flux_s, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:2077
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1394
virtual ConvergenceReport getWellConvergence(const Simulator &simulator, const WellStateType &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:1170
void checkConvergenceExtraEqs(const std::vector< Scalar > &res, ConvergenceReport &report) const
Definition: StandardWell_impl.hpp:2125
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:1670
bool openCrossFlowAvoidSingularity(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1108
std::optional< Scalar > computeBhpAtThpLimitInj(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2283
bool allDrawDownWrongDirection(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1066
EvalWell pskin(const Scalar throughput, const EvalWell &water_velocity, const EvalWell &poly_inj_conc, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1949
static constexpr int numWellConservationEq
Definition: StandardWell.hpp:95
int setPrimaryVars(typename std::vector< Scalar >::const_iterator it) override
Definition: StandardWell_impl.hpp:2535
void checkOperabilityUnderTHPLimit(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1017
void computeWellPotentials(const Simulator &simulator, const WellStateType &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) override
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1716
bool computeWellPotentialsImplicit(const Simulator &ebos_simulator, const WellStateType &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1591
void updateWaterThroughput(const double dt, WellStateType &well_state) const override
Definition: StandardWell_impl.hpp:2023
void checkOperabilityUnderBHPLimit(const WellStateType &well_state, const Simulator &simulator, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:948
void recoverWellSolutionAndUpdateWellState(const Simulator &simulator, const BVector &x, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1426
void assembleWellEqWithoutIterationImpl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:367
bool iterateWellEqWithControl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:2315
EvalWell pskinwater(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1917
void solveEqAndUpdateWellState(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1358
void handleInjectivityRate(const Simulator &simulator, const int perf, std::vector< EvalWell > &cq_s) const
Definition: StandardWell_impl.hpp:2054
virtual void init(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
std::vector< Scalar > computeWellPotentialWithTHP(const Simulator &ebosSimulator, DeferredLogger &deferred_logger, const WellStateType &well_state) const
Definition: StandardWell_impl.hpp:1556
void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator< Scalar > &wellPICalc, WellStateType &well_state, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1214
void updatePrimaryVariables(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1791
Scalar getRefDensity() const override
Definition: StandardWell_impl.hpp:1814
Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const override
Definition: StandardWell_impl.hpp:1776
EvalWell getQs(const int compIdx) const
Returns scaled rate for a component.
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
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.
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.
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
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.
Definition: WellConvergence.hpp:38
const int num_conservation_quantities_
Definition: WellInterfaceGeneric.hpp:310
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:300
void resetDampening()
Definition: WellInterfaceGeneric.hpp:243
std::pair< bool, bool > computeWellPotentials(std::vector< Scalar > &well_potentials, const WellStateType &well_state)
void prepareForPotentialCalculations(const SummaryState &summary_state, WellStateType &well_state, Well::InjectionControls &inj_controls, Well::ProductionControls &prod_controls) const
Definition: WellInterfaceIndices.hpp:34
Definition: WellInterface.hpp:76
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:81
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:103
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:1899
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:86
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:82
Dune::BlockVector< VectorBlockType > BVector
Definition: WellInterface.hpp:96
typename Base::ModelParameters ModelParameters
Definition: WellInterface.hpp:109
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:83
bool solveWellWithOperabilityCheck(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:563
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:66
const SingleWellState< Scalar, IndexTraits > & well(std::size_t well_index) const
Definition: WellState.hpp:289
std::vector< Scalar > & wellRates(std::size_t well_index)
One rate per well and phase.
Definition: WellState.hpp:254
@ NONE
Definition: DeferredLogger.hpp:46
Definition: blackoilboundaryratevector.hh:39
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