23#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
24#define OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
29#include <opm/common/Exceptions.hpp>
31#include <opm/input/eclipse/Units/Units.hpp>
33#include <opm/material/densead/EvaluationFormat.hpp>
41#include <fmt/format.h>
50template<
class dValue,
class Value>
51auto dValueError(
const dValue& d,
52 const std::string& name,
53 const std::string& methodName,
56 const Value& pressure)
58 return fmt::format(
"Problematic d value {} obtained for well {}"
59 " during {} calculations with rs {}"
60 ", rv {} and pressure {}."
61 " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
62 " for this connection.", d, name, methodName, Rs, Rv, pressure);
70 template<
typename TypeTag>
77 const int pvtRegionIdx,
78 const int num_components,
80 const int index_of_well,
81 const std::vector<PerforationData>& perf_data)
82 :
Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
93 template<
typename TypeTag>
97 const std::vector<double>& depth_arg,
98 const double gravity_arg,
100 const std::vector< Scalar >& B_avg,
101 const bool changed_to_open_this_step)
103 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
104 this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
111 template<
typename TypeTag>
115 this->primary_variables_.init();
122 template<
typename TypeTag>
123 template<
class Value>
127 const std::vector<Value>& mob,
129 const std::vector<Scalar>& Tw,
132 std::vector<Value>& cq_s,
136 auto obtain = [
this](
const Eval& value)
138 if constexpr (std::is_same_v<Value, Scalar>) {
139 static_cast<void>(
this);
140 return getValue(value);
142 return this->extendEval(value);
145 auto obtainN = [](
const auto& value)
147 if constexpr (std::is_same_v<Value, Scalar>) {
148 return getValue(value);
153 auto zeroElem = [
this]()
155 if constexpr (std::is_same_v<Value, Scalar>) {
156 static_cast<void>(
this);
159 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
163 const auto& fs = intQuants.fluidState();
164 const Value pressure = obtain(this->getPerfCellPressure(fs));
165 const Value rs = obtain(fs.Rs());
166 const Value rv = obtain(fs.Rv());
167 const Value rvw = obtain(fs.Rvw());
168 const Value rsw = obtain(fs.Rsw());
170 std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
171 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
172 if (!FluidSystem::phaseIsActive(phaseIdx)) {
175 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
176 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
178 if constexpr (has_solvent) {
179 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
182 if constexpr (has_zFraction) {
183 if (this->isInjector()) {
184 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
185 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
186 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
190 Value skin_pressure = zeroElem();
192 if (this->isInjector()) {
193 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
194 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
199 std::vector<Value> cmix_s(this->numComponents(), zeroElem());
200 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
201 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
222 template<
typename TypeTag>
223 template<
class Value>
227 const Value& pressure,
233 std::vector<Value>& b_perfcells_dense,
234 const std::vector<Scalar>& Tw,
237 const Value& skin_pressure,
238 const std::vector<Value>& cmix_s,
239 std::vector<Value>& cq_s,
244 const Value well_pressure =
bhp + this->connections_.pressure_diff(perf);
245 Value drawdown = pressure - well_pressure;
246 if (this->isInjector()) {
247 drawdown += skin_pressure;
253 if (!allow_cf && this->isInjector()) {
258 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
259 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
260 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
263 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
264 gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw);
265 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
266 gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw);
270 if (!allow_cf && this->isProducer()) {
275 Value total_mob_dense = mob[0];
276 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
277 total_mob_dense += mob[componentIdx];
281 Value volumeRatio =
bhp * 0.0;
283 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
284 disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
285 cmix_s, b_perfcells_dense, deferred_logger);
287 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
288 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
289 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
293 if constexpr (Indices::enableSolvent) {
294 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
297 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
298 gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
299 cmix_s, b_perfcells_dense, deferred_logger);
302 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
303 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
304 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
306 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
307 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
308 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
313 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
314 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
315 Value cqt_is = cqt_i / volumeRatio;
316 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
320 if (this->isProducer()) {
321 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
322 gasOilPerfRateInj(cq_s, perf_rates,
323 rv, rs, pressure, rvw, deferred_logger);
325 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
327 gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
328 pressure, deferred_logger);
335 template<
typename TypeTag>
340 const Well::InjectionControls& inj_controls,
341 const Well::ProductionControls& prod_controls,
348 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
351 this->linSys_.clear();
353 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
354 prod_controls, well_state,
355 group_state, deferred_logger);
361 template<
typename TypeTag>
366 const Well::InjectionControls& inj_controls,
367 const Well::ProductionControls& prod_controls,
373 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
374 const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
376 auto& ws = well_state.
well(this->index_of_well_);
377 ws.phase_mixing_rates.fill(0.0);
380 const int np = this->number_of_phases_;
382 std::vector<RateVector> connectionRates = this->connectionRates_;
384 auto& perf_data = ws.perf_data;
385 auto& perf_rates = perf_data.phase_rates;
386 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
388 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
389 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
390 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
391 calculateSinglePerf(simulator, perf, well_state, connectionRates,
392 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
395 if constexpr (has_polymer && Base::has_polymermw) {
396 if (this->isInjector()) {
397 handleInjectivityEquations(simulator, well_state, perf,
398 water_flux_s, deferred_logger);
401 const int cell_idx = this->well_cells_[perf];
402 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
404 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
406 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
409 assemblePerforationEq(cq_s_effective,
412 this->primary_variables_.numWellEq(),
416 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
417 auto& perf_rate_solvent = perf_data.solvent_rates;
418 perf_rate_solvent[perf] = cq_s[componentIdx].value();
420 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
424 if constexpr (has_zFraction) {
426 assembleZFracEq(cq_s_zfrac_effective,
428 this->primary_variables_.numWellEq(),
433 this->connectionRates_ = connectionRates;
438 const auto& comm = this->parallel_well_info_.communication();
439 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
443 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
446 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
449 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
450 if (FluidSystem::numActivePhases() > 1) {
452 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
453 this->F0_[componentIdx]) * volume / dt;
455 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
457 assembleSourceEq(resWell_loc,
459 this->primary_variables_.numWellEq(),
463 const auto& summaryState = simulator.vanguard().summaryState();
464 const Schedule& schedule = simulator.vanguard().schedule();
466 assembleControlEq(well_state, group_state,
467 schedule, summaryState,
468 inj_controls, prod_controls,
469 this->primary_variables_,
470 this->connections_.rho(),
477 this->linSys_.invert();
479 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
486 template<
typename TypeTag>
492 std::vector<RateVector>& connectionRates,
493 std::vector<EvalWell>& cq_s,
498 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
499 const EvalWell&
bhp = this->primary_variables_.eval(Bhp);
500 const int cell_idx = this->well_cells_[perf];
501 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
502 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
503 getMobility(simulator, perf, mob, deferred_logger);
506 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
507 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
508 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
509 computePerfRate(intQuants, mob,
bhp, Tw, perf, allow_cf,
510 cq_s, perf_rates, deferred_logger);
512 auto& ws = well_state.
well(this->index_of_well_);
513 auto& perf_data = ws.perf_data;
514 if constexpr (has_polymer && Base::has_polymermw) {
515 if (this->isInjector()) {
518 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
519 water_flux_s = cq_s[water_comp_idx];
522 handleInjectivityRate(simulator, perf, cq_s);
527 if (this->isProducer()) {
528 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.
dis_gas;
529 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.
dis_gas_in_water;
530 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.
vap_oil;
531 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.
vap_wat;
532 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.
dis_gas;
533 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.
dis_gas_in_water;
534 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.
vap_oil;
535 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.
vap_wat;
538 if constexpr (has_energy) {
539 connectionRates[perf][Indices::contiEnergyEqIdx] =
540 connectionRateEnergy(simulator.problem().maxOilSaturation(cell_idx),
541 cq_s, intQuants, deferred_logger);
544 if constexpr (has_polymer) {
545 std::variant<Scalar,EvalWell> polymerConcentration;
546 if (this->isInjector()) {
547 polymerConcentration = this->wpolymer();
549 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
550 intQuants.polymerViscosityCorrection());
553 [[maybe_unused]]
EvalWell cq_s_poly;
554 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
556 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
557 cq_s, polymerConcentration);
559 if constexpr (Base::has_polymermw) {
560 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
561 perf, connectionRates, deferred_logger);
565 if constexpr (has_foam) {
566 std::variant<Scalar,EvalWell> foamConcentration;
567 if (this->isInjector()) {
568 foamConcentration = this->wfoam();
570 foamConcentration = this->extendEval(intQuants.foamConcentration());
572 connectionRates[perf][Indices::contiFoamEqIdx] =
573 this->connections_.connectionRateFoam(cq_s, foamConcentration,
574 FoamModule::transportPhase(),
578 if constexpr (has_zFraction) {
579 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
580 if (this->isInjector()) {
581 solventConcentration = this->wsolvent();
583 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
584 this->extendEval(intQuants.yVolume())};
586 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
587 cq_s_zfrac_effective) =
588 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
590 solventConcentration);
593 if constexpr (has_brine) {
594 std::variant<Scalar,EvalWell> saltConcentration;
595 if (this->isInjector()) {
596 saltConcentration = this->wsalt();
598 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
601 connectionRates[perf][Indices::contiBrineEqIdx] =
602 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
607 if constexpr (has_micp) {
608 std::variant<Scalar,EvalWell> microbialConcentration;
609 std::variant<Scalar,EvalWell> oxygenConcentration;
610 std::variant<Scalar,EvalWell> ureaConcentration;
611 if (this->isInjector()) {
612 microbialConcentration = this->wmicrobes();
613 oxygenConcentration = this->woxygen();
614 ureaConcentration = this->wurea();
616 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
617 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
618 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
620 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
621 connectionRates[perf][Indices::contiOxygenEqIdx],
622 connectionRates[perf][Indices::contiUreaEqIdx]) =
623 this->connections_.connectionRatesMICP(cq_s,
624 microbialConcentration,
630 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
635 template<
typename TypeTag>
636 template<
class Value>
641 std::vector<Value>& mob,
644 auto obtain = [
this](
const Eval& value)
646 if constexpr (std::is_same_v<Value, Scalar>) {
647 static_cast<void>(
this);
648 return getValue(value);
650 return this->extendEval(value);
654 obtain, deferred_logger);
657 if constexpr (has_polymer) {
658 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
659 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
664 if constexpr (!Base::has_polymermw) {
665 if constexpr (std::is_same_v<Value, Scalar>) {
666 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
667 for (std::size_t i = 0; i < mob.size(); ++i) {
668 mob_eval[i].setValue(mob[i]);
670 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
671 for (std::size_t i = 0; i < mob.size(); ++i) {
672 mob[i] = getValue(mob_eval[i]);
675 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
682 const double bhp = this->primary_variables_.
value(Bhp);
683 const double perf_press =
bhp + this->connections_.pressure_diff(perf);
684 const double multiplier = this->getInjMult(perf,
bhp, perf_press);
685 for (std::size_t i = 0; i < mob.size(); ++i) {
686 mob[i] *= multiplier;
692 template<
typename TypeTag>
700 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
702 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
703 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
705 updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
706 Base::calculateReservoirRates(well_state.
well(this->index_of_well_));
713 template<
typename TypeTag>
717 const bool stop_or_zero_rate_target,
720 const double dFLimit = this->param_.dwell_fraction_max_;
721 const double dBHPLimit = this->param_.dbhp_max_rel_;
722 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
725 if constexpr (Base::has_polymermw) {
726 this->primary_variables_.updateNewtonPolyMW(dwells);
729 this->primary_variables_.checkFinite(deferred_logger);
736 template<
typename TypeTag>
741 const SummaryState& summary_state,
744 this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
747 if constexpr (Base::has_polymermw) {
748 this->primary_variables_.copyToWellStatePolyMW(well_state);
756 template<
typename TypeTag>
764 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
765 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
767 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
768 std::vector<Scalar> mob(this->num_components_, 0.0);
769 getMobility(simulator, perf, mob, deferred_logger);
771 const int cell_idx = this->well_cells_[perf];
772 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
773 const auto& fs = int_quantities.fluidState();
775 double p_r = this->getPerfCellPressure(fs).value();
778 std::vector<double> b_perf(this->num_components_);
779 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
780 if (!FluidSystem::phaseIsActive(phase)) {
783 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
784 b_perf[comp_idx] = fs.invB(phase).value();
786 if constexpr (has_solvent) {
787 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
791 const double h_perf = this->connections_.pressure_diff(perf);
792 const double pressure_diff = p_r - h_perf;
797 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
798 deferred_logger.
debug(
"CROSSFLOW_IPR",
799 "cross flow found when updateIPR for well " + name()
800 +
" . The connection is ignored in IPR calculations");
806 double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
807 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
808 const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
809 std::vector<double> ipr_a_perf(this->ipr_a_.size());
810 std::vector<double> ipr_b_perf(this->ipr_b_.size());
811 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
812 const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
813 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
814 ipr_b_perf[comp_idx] += tw_mob;
818 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
819 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
820 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
821 const double rs = (fs.Rs()).value();
822 const double rv = (fs.Rv()).value();
824 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
825 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
827 ipr_a_perf[gas_comp_idx] += dis_gas_a;
828 ipr_a_perf[oil_comp_idx] += vap_oil_a;
830 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
831 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
833 ipr_b_perf[gas_comp_idx] += dis_gas_b;
834 ipr_b_perf[oil_comp_idx] += vap_oil_b;
837 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
838 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
839 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
842 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
843 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
846 template<
typename TypeTag>
860 auto rates = well_state.
well(this->index_of_well_).surface_rates;
862 for (std::size_t p = 0; p < rates.size(); ++p) {
863 zero_rates &= rates[p] == 0.0;
865 auto& ws = well_state.
well(this->index_of_well_);
867 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
868 deferred_logger.
debug(msg);
880 const auto& group_state = simulator.problem().wellModel().groupState();
882 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
883 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
885 auto inj_controls = Well::InjectionControls(0);
886 auto prod_controls = Well::ProductionControls(0);
887 prod_controls.addControl(Well::ProducerCMode::BHP);
888 prod_controls.bhp_limit = well_state.
well(this->index_of_well_).bhp;
891 const auto cmode = ws.production_cmode;
892 ws.production_cmode = Well::ProducerCMode::BHP;
893 const double dt = simulator.timeStepSize();
894 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
896 const size_t nEq = this->primary_variables_.numWellEq();
900 for (
size_t i=0; i < nEq; ++i){
906 x_well[0].resize(nEq);
907 this->linSys_.solve(rhs, x_well);
909 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
910 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
911 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
912 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
914 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
916 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
919 ws.production_cmode = cmode;
922 template<
typename TypeTag>
929 const auto& summaryState = simulator.vanguard().summaryState();
933 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
934 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
937 double total_ipr_mass_rate = 0.0;
938 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
940 if (!FluidSystem::phaseIsActive(phaseIdx)) {
944 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
945 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
947 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
948 total_ipr_mass_rate += ipr_rate * rho;
950 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
951 this->operability_status_.operable_under_only_bhp_limit =
false;
955 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
959 std::vector<double> well_rates_bhp_limit;
960 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
962 this->adaptRatesForVFP(well_rates_bhp_limit);
963 const double thp_limit = this->getTHPConstraint(summaryState);
966 this->connections_.rho(),
967 this->getALQ(well_state),
970 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
971 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
982 this->operability_status_.operable_under_only_bhp_limit =
true;
983 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
991 template<
typename TypeTag>
998 const auto& summaryState = simulator.vanguard().summaryState();
999 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1000 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1003 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1006 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1007 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1009 const double thp_limit = this->getTHPConstraint(summaryState);
1010 if (this->isProducer() && *obtain_bhp < thp_limit) {
1011 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1012 +
" bars is SMALLER than thp limit "
1014 +
" bars as a producer for well " + name();
1015 deferred_logger.
debug(msg);
1017 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1018 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1019 +
" bars is LARGER than thp limit "
1021 +
" bars as a injector for well " + name();
1022 deferred_logger.
debug(msg);
1025 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1026 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1027 if (!this->wellIsStopped()) {
1028 const double thp_limit = this->getTHPConstraint(summaryState);
1029 deferred_logger.
debug(
" could not find bhp value at thp limit "
1031 +
" bar for well " + name() +
", the well might need to be closed ");
1040 template<
typename TypeTag>
1045 bool all_drawdown_wrong_direction =
true;
1047 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1048 const int cell_idx = this->well_cells_[perf];
1049 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1050 const auto& fs = intQuants.fluidState();
1052 const double pressure = this->getPerfCellPressure(fs).value();
1053 const double bhp = this->primary_variables_.eval(Bhp).
value();
1056 const double well_pressure =
bhp + this->connections_.pressure_diff(perf);
1057 const double drawdown = pressure - well_pressure;
1062 if ( (drawdown < 0. && this->isInjector()) ||
1063 (drawdown > 0. && this->isProducer()) ) {
1064 all_drawdown_wrong_direction =
false;
1069 const auto& comm = this->parallel_well_info_.communication();
1070 if (comm.size() > 1)
1072 all_drawdown_wrong_direction =
1073 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1076 return all_drawdown_wrong_direction;
1082 template<
typename TypeTag>
1089 const double bhp = well_state.
well(this->index_of_well_).bhp;
1090 std::vector<double> well_rates;
1091 computeWellRatesWithBhp(simulator,
bhp, well_rates, deferred_logger);
1093 const double sign = (this->isProducer()) ? -1. : 1.;
1094 const double threshold = sign * std::numeric_limits<double>::min();
1096 bool can_produce_inject =
false;
1097 for (
const auto value : well_rates) {
1098 if (this->isProducer() && value < threshold) {
1099 can_produce_inject =
true;
1101 }
else if (this->isInjector() && value > threshold) {
1102 can_produce_inject =
true;
1107 if (!can_produce_inject) {
1108 deferred_logger.
debug(
" well " + name() +
" CANNOT produce or inejct ");
1111 return can_produce_inject;
1118 template<
typename TypeTag>
1123 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1129 template<
typename TypeTag>
1136 std::function<
Scalar(
int,
int)> getTemperature =
1137 [&simulator](
int cell_idx,
int phase_idx)
1139 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().temperature(phase_idx).value();
1141 std::function<
Scalar(
int)> getSaltConcentration =
1142 [&simulator](
int cell_idx)
1144 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().saltConcentration().value();
1146 std::function<int(
int)> getPvtRegionIdx =
1147 [&simulator](
int cell_idx)
1149 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().pvtRegionIndex();
1151 std::function<
Scalar(
int)> getInvFac =
1152 [&simulator](
int cell_idx)
1154 return simulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1156 std::function<
Scalar(
int)> getSolventDensity =
1157 [&simulator](
int cell_idx)
1159 return simulator.model().intensiveQuantities(cell_idx, 0).solventRefDensity();
1162 this->connections_.computePropertiesForPressures(well_state,
1164 getSaltConcentration,
1175 template<
typename TypeTag>
1180 const std::vector<double>& B_avg,
1182 const bool relax_tolerance)
const
1186 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1188 double tol_wells = this->param_.tolerance_wells_;
1190 constexpr double stopped_factor = 1.e-4;
1192 constexpr double dynamic_thp_factor = 1.e-1;
1193 if (this->stopppedOrZeroRateTarget(summary_state, well_state)) {
1194 tol_wells = tol_wells*stopped_factor;
1195 }
else if (this->getDynamicThpLimit()) {
1196 tol_wells = tol_wells*dynamic_thp_factor;
1199 std::vector<double> res;
1202 this->param_.max_residual_allowed_,
1204 this->param_.relaxed_tolerance_flow_well_,
1206 this->wellIsStopped(),
1210 checkConvergenceExtraEqs(res, report);
1219 template<
typename TypeTag>
1227 auto fluidState = [&simulator,
this](
const int perf)
1229 const auto cell_idx = this->well_cells_[perf];
1230 return simulator.model()
1231 .intensiveQuantities(cell_idx, 0).fluidState();
1234 const int np = this->number_of_phases_;
1235 auto setToZero = [np](
double* x) ->
void
1237 std::fill_n(x, np, 0.0);
1240 auto addVector = [np](
const double* src,
double* dest) ->
void
1242 std::transform(src, src + np, dest, dest, std::plus<>{});
1245 auto& ws = well_state.
well(this->index_of_well_);
1246 auto& perf_data = ws.perf_data;
1247 auto* wellPI = ws.productivity_index.data();
1248 auto* connPI = perf_data.prod_index.data();
1252 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1253 auto subsetPerfID = 0;
1255 for (
const auto& perf : *this->perf_data_) {
1256 auto allPerfID = perf.ecl_index;
1258 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1263 std::vector<Scalar> mob(this->num_components_, 0.0);
1264 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1266 const auto& fs = fluidState(subsetPerfID);
1269 if (this->isInjector()) {
1270 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1271 mob, connPI, deferred_logger);
1274 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1277 addVector(connPI, wellPI);
1284 const auto& comm = this->parallel_well_info_.communication();
1285 if (comm.size() > 1) {
1286 comm.sum(wellPI, np);
1289 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1290 "Internal logic error in processing connections for PI/II");
1295 template<
typename TypeTag>
1303 std::function<
Scalar(
int,
int)> invB =
1304 [&simulator](
int cell_idx,
int phase_idx)
1306 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().invB(phase_idx).value();
1308 std::function<
Scalar(
int,
int)> mobility =
1309 [&simulator](
int cell_idx,
int phase_idx)
1311 return simulator.model().intensiveQuantities(cell_idx, 0).mobility(phase_idx).value();
1313 std::function<
Scalar(
int)> invFac =
1314 [&simulator](
int cell_idx)
1316 return simulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1318 std::function<
Scalar(
int)> solventMobility =
1319 [&simulator](
int cell_idx)
1321 return simulator.model().intensiveQuantities(cell_idx, 0).solventMobility().value();
1324 this->connections_.computeProperties(well_state,
1337 template<
typename TypeTag>
1348 computePropertiesForWellConnectionPressures(simulator, well_state, props);
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(summary_state, dx_well, well_state, deferred_logger);
1379 template<
typename TypeTag>
1386 const auto& summary_state = simulator.vanguard().summaryState();
1387 updatePrimaryVariables(summary_state, well_state, deferred_logger);
1388 initPrimaryVariablesEvaluation();
1389 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1390 this->computeAccumWell();
1395 template<
typename TypeTag>
1400 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1402 if (this->param_.matrix_add_well_contributions_)
1408 this->linSys_.apply(x, Ax);
1414 template<
typename TypeTag>
1419 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1421 this->linSys_.apply(r);
1427 template<
typename TypeTag>
1435 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1438 xw[0].resize(this->primary_variables_.numWellEq());
1440 this->linSys_.recoverSolutionWell(x, xw);
1441 updateWellState(summary_state, xw, well_state, deferred_logger);
1447 template<
typename TypeTag>
1452 std::vector<double>& well_flux,
1456 const int np = this->number_of_phases_;
1457 well_flux.resize(np, 0.0);
1459 const bool allow_cf = this->getAllowCrossFlow();
1461 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1462 const int cell_idx = this->well_cells_[perf];
1463 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1465 std::vector<Scalar> mob(this->num_components_, 0.);
1466 getMobility(simulator, perf, mob, deferred_logger);
1467 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
1468 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1469 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1471 std::vector<Scalar> cq_s(this->num_components_, 0.);
1473 computePerfRate(intQuants, mob,
bhp, Tw, perf, allow_cf,
1474 cq_s, perf_rates, deferred_logger);
1476 for(
int p = 0; p < np; ++p) {
1477 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1481 if constexpr (has_solvent) {
1483 assert(pu.phase_used[Gas]);
1484 const int gas_pos = pu.phase_pos[Gas];
1485 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1488 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1493 template<
typename TypeTag>
1498 std::vector<double>& well_flux,
1508 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1509 const auto& group_state = simulator.problem().wellModel().groupState();
1512 const auto& summary_state = simulator.vanguard().summaryState();
1513 auto inj_controls = well_copy.
well_ecl_.isInjector()
1514 ? well_copy.
well_ecl_.injectionControls(summary_state)
1515 : Well::InjectionControls(0);
1516 auto prod_controls = well_copy.
well_ecl_.isProducer()
1517 ? well_copy.
well_ecl_.productionControls(summary_state) :
1518 Well::ProductionControls(0);
1521 auto& ws = well_state_copy.
well(this->index_of_well_);
1523 inj_controls.bhp_limit =
bhp;
1524 ws.injection_cmode = Well::InjectorCMode::BHP;
1526 prod_controls.bhp_limit =
bhp;
1527 ws.production_cmode = Well::ProducerCMode::BHP;
1532 const int np = this->number_of_phases_;
1533 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1534 for (
int phase = 0; phase < np; ++phase){
1535 well_state_copy.
wellRates(this->index_of_well_)[phase]
1536 = sign * ws.well_potentials[phase];
1542 const double dt = simulator.timeStepSize();
1543 const bool converged = well_copy.
iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1545 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1546 " potentials are computed based on unconverged solution";
1547 deferred_logger.
debug(msg);
1558 template<
typename TypeTag>
1565 std::vector<double> potentials(this->number_of_phases_, 0.0);
1566 const auto& summary_state = simulator.vanguard().summaryState();
1568 const auto& well = this->well_ecl_;
1569 if (well.isInjector()){
1570 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1571 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1572 if (bhp_at_thp_limit) {
1573 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1574 computeWellRatesWithBhp(simulator,
bhp, potentials, deferred_logger);
1576 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1577 "Failed in getting converged thp based potential calculation for well "
1578 + name() +
". Instead the bhp based value is used");
1579 const double bhp = controls.bhp_limit;
1580 computeWellRatesWithBhp(simulator,
bhp, potentials, deferred_logger);
1583 computeWellRatesWithThpAlqProd(
1584 simulator, summary_state,
1585 deferred_logger, potentials, this->getALQ(well_state)
1592 template<
typename TypeTag>
1596 std::vector<double>& well_potentials,
1605 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1606 const auto& group_state = simulator.problem().wellModel().groupState();
1607 auto& ws = well_state_copy.
well(this->index_of_well_);
1610 const auto& summary_state = simulator.vanguard().summaryState();
1611 auto inj_controls = well_copy.
well_ecl_.isInjector()
1612 ? well_copy.
well_ecl_.injectionControls(summary_state)
1613 : Well::InjectionControls(0);
1614 auto prod_controls = well_copy.
well_ecl_.isProducer()
1615 ? well_copy.
well_ecl_.productionControls(summary_state) :
1616 Well::ProductionControls(0);
1622 const int np = this->number_of_phases_;
1623 bool trivial =
true;
1624 for (
int phase = 0; phase < np; ++phase){
1625 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1628 const double sign = well_copy.
well_ecl_.isInjector() ? 1.0 : -1.0;
1629 for (
int phase = 0; phase < np; ++phase) {
1630 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1635 const double dt = simulator.timeStepSize();
1637 bool converged =
false;
1638 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1639 converged = well_copy.
solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1641 converged = well_copy.
iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1645 well_potentials.clear();
1646 well_potentials.resize(np, 0.0);
1647 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
1649 well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
1655 template<
typename TypeTag>
1659 const SummaryState &summary_state,
1661 std::vector<double> &potentials,
1665 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1666 simulator, summary_state, alq, deferred_logger);
1667 if (bhp_at_thp_limit) {
1668 const auto& controls = this->well_ecl_.productionControls(summary_state);
1669 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1670 computeWellRatesWithBhp(simulator,
bhp, potentials, deferred_logger);
1673 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1674 "Failed in getting converged thp based potential calculation for well "
1675 + name() +
". Instead the bhp based value is used");
1676 const auto& controls = this->well_ecl_.productionControls(summary_state);
1677 bhp = controls.bhp_limit;
1678 computeWellRatesWithBhp(simulator,
bhp, potentials, deferred_logger);
1683 template<
typename TypeTag>
1687 const SummaryState& summary_state,
1689 std::vector<double>& potentials,
1693 computeWellRatesAndBhpWithThpAlqProd(simulator,
1700 template<
typename TypeTag>
1705 std::vector<double>& well_potentials,
1708 const auto [compute_potential, bhp_controlled_well] =
1711 if (!compute_potential) {
1715 bool converged_implicit =
false;
1716 if (this->param_.local_well_solver_control_switching_) {
1717 converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
1719 if (!converged_implicit) {
1721 const auto& summaryState = simulator.vanguard().summaryState();
1722 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1732 const auto& ws = well_state.
well(this->index_of_well_);
1733 if (this->isInjector())
1734 bhp = std::max(ws.bhp,
bhp);
1736 bhp = std::min(ws.bhp,
bhp);
1738 assert(std::abs(
bhp) != std::numeric_limits<double>::max());
1739 computeWellRatesWithBhpIterations(simulator,
bhp, well_potentials, deferred_logger);
1742 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1746 this->checkNegativeWellPotentials(well_potentials,
1747 this->param_.check_well_operability_,
1757 template<
typename TypeTag>
1761 const int openConnIdx)
const
1763 return (openConnIdx < 0)
1765 : this->connections_.rho(openConnIdx);
1772 template<
typename TypeTag>
1779 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1781 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
1782 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1785 if constexpr (Base::has_polymermw) {
1786 this->primary_variables_.updatePolyMW(well_state);
1789 this->primary_variables_.checkFinite(deferred_logger);
1795 template<
typename TypeTag>
1800 return this->connections_.rho();
1806 template<
typename TypeTag>
1811 std::vector<EvalWell>& mob,
1814 const int cell_idx = this->well_cells_[perf];
1815 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1816 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1820 if (this->isInjector()) {
1822 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1823 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1824 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1827 if (PolymerModule::hasPlyshlog()) {
1830 if (this->isInjector() && this->wpolymer() == 0.) {
1835 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1836 const EvalWell&
bhp = this->primary_variables_.eval(Bhp);
1838 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1840 double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quant, cell_idx);
1841 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1842 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1843 computePerfRate(int_quant, mob,
bhp, Tw, perf, allow_cf, cq_s,
1844 perf_rates, deferred_logger);
1846 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1847 const auto& material_law_manager = simulator.problem().materialLawManager();
1848 const auto& scaled_drainage_info =
1849 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1850 const double swcr = scaled_drainage_info.Swcr;
1851 const EvalWell poro = this->extendEval(int_quant.porosity());
1852 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1854 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1855 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1856 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1858 if (PolymerModule::hasShrate()) {
1861 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1863 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1864 int_quant.pvtRegionIndex(),
1867 mob[waterCompIdx] /= shear_factor;
1871 template<
typename TypeTag>
1875 this->linSys_.extract(jacobian);
1879 template <
typename TypeTag>
1883 const int pressureVarIndex,
1884 const bool use_well_weights,
1887 this->linSys_.extractCPRPressureMatrix(jacobian,
1898 template<
typename TypeTag>
1905 if constexpr (Base::has_polymermw) {
1906 const int water_table_id = this->polymerWaterTable_();
1907 if (water_table_id <= 0) {
1909 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1912 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1913 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1915 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1916 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1920 fmt::format(
"Polymermw is not activated, while injecting "
1921 "skin pressure is requested for well {}", name()),
1930 template<
typename TypeTag>
1933 pskin(
const double throughput,
1938 if constexpr (Base::has_polymermw) {
1939 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
1940 const EvalWell water_velocity_abs = abs(water_velocity);
1941 if (poly_inj_conc == 0.) {
1942 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1944 const int polymer_table_id = this->polymerTable_();
1945 if (polymer_table_id <= 0) {
1947 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1950 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1951 const double reference_concentration = skprpolytable.refConcentration;
1952 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1954 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1955 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1956 if (poly_inj_conc == reference_concentration) {
1957 return sign * pskin_poly;
1960 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1961 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1962 return sign * pskin;
1965 fmt::format(
"Polymermw is not activated, while injecting "
1966 "skin pressure is requested for well {}", name()),
1975 template<
typename TypeTag>
1982 if constexpr (Base::has_polymermw) {
1983 const int table_id = this->polymerInjTable_();
1984 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
1985 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1986 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
1987 if (this->wpolymer() == 0.) {
1988 return molecular_weight;
1990 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
1991 return molecular_weight;
1994 fmt::format(
"Polymermw is not activated, while injecting "
1995 "polymer molecular weight is requested for well {}", name()),
2004 template<
typename TypeTag>
2009 if constexpr (Base::has_polymermw) {
2010 if (this->isInjector()) {
2011 auto& ws = well_state.
well(this->index_of_well_);
2012 auto& perf_water_throughput = ws.perf_data.water_throughput;
2013 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2014 const double perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
2016 if (perf_water_vel > 0.) {
2017 perf_water_throughput[perf] += perf_water_vel * dt;
2028 template<
typename TypeTag>
2033 std::vector<EvalWell>& cq_s)
const
2035 const int cell_idx = this->well_cells_[perf];
2036 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2037 const auto& fs = int_quants.fluidState();
2038 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2039 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2040 const int wat_vel_index = Bhp + 1 + perf;
2041 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2045 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2051 template<
typename TypeTag>
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 EvalWell water_flux_r = water_flux_s / b_w;
2065 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2066 const EvalWell water_velocity = water_flux_r / area;
2067 const int wat_vel_index = Bhp + 1 + perf;
2070 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2072 const auto& ws = well_state.
well(this->index_of_well_);
2073 const auto& perf_data = ws.perf_data;
2074 const auto& perf_water_throughput = perf_data.water_throughput;
2075 const double throughput = perf_water_throughput[perf];
2076 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2078 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2079 poly_conc.setValue(this->wpolymer());
2082 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2083 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2086 assembleInjectivityEq(eq_pskin,
2091 this->primary_variables_.numWellEq(),
2099 template<
typename TypeTag>
2108 if constexpr (Base::has_polymermw) {
2110 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2118 template<
typename TypeTag>
2125 std::vector<RateVector>& connectionRates,
2130 if (this->isInjector()) {
2131 const int wat_vel_index = Bhp + 1 + perf;
2132 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2133 if (water_velocity > 0.) {
2134 const auto& ws = well_state.
well(this->index_of_well_);
2135 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2136 const double throughput = perf_water_throughput[perf];
2137 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2138 cq_s_polymw *= molecular_weight;
2144 }
else if (this->isProducer()) {
2145 if (cq_s_polymw < 0.) {
2146 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2153 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2161 template<
typename TypeTag>
2162 std::optional<double>
2166 const SummaryState& summary_state,
2169 return computeBhpAtThpLimitProdWithAlq(simulator,
2171 this->getALQ(well_state),
2175 template<
typename TypeTag>
2176 std::optional<double>
2179 const SummaryState& summary_state,
2180 const double alq_value,
2184 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2190 std::vector<double> rates(3);
2191 computeWellRatesWithBhp(simulator,
bhp, rates, deferred_logger);
2192 this->adaptRatesForVFP(rates);
2196 double max_pressure = 0.0;
2197 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2198 const int cell_idx = this->well_cells_[perf];
2199 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2200 const auto& fs = int_quants.fluidState();
2201 double pressure_cell = this->getPerfCellPressure(fs).value();
2202 max_pressure = std::max(max_pressure, pressure_cell);
2207 this->connections_.rho(),
2209 this->getTHPConstraint(summary_state),
2213 auto v = frates(*bhpAtLimit);
2214 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2220 auto fratesIter = [
this, &simulator, &deferred_logger](
const double bhp) {
2224 std::vector<double> rates(3);
2225 computeWellRatesWithBhpIterations(simulator,
bhp, rates, deferred_logger);
2226 this->adaptRatesForVFP(rates);
2233 this->connections_.rho(),
2235 this->getTHPConstraint(summary_state),
2241 auto v = frates(*bhpAtLimit);
2242 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2248 return std::nullopt;
2253 template<
typename TypeTag>
2254 std::optional<double>
2257 const SummaryState& summary_state,
2261 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2267 std::vector<double> rates(3);
2268 computeWellRatesWithBhp(simulator,
bhp, rates, deferred_logger);
2274 this->connections_.rho(),
2285 template<
typename TypeTag>
2290 const Well::InjectionControls& inj_controls,
2291 const Well::ProductionControls& prod_controls,
2296 const int max_iter = this->param_.max_inner_iter_wells_;
2299 bool relax_convergence =
false;
2300 this->regularize_ =
false;
2301 const auto& summary_state = simulator.vanguard().summaryState();
2303 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2305 if (it > this->param_.strict_inner_iter_wells_) {
2306 relax_convergence =
true;
2307 this->regularize_ =
true;
2310 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2312 converged = report.converged();
2318 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2325 initPrimaryVariablesEvaluation();
2326 }
while (it < max_iter);
2332 template<
typename TypeTag>
2337 const Well::InjectionControls& inj_controls,
2338 const Well::ProductionControls& prod_controls,
2342 const bool fixed_control ,
2343 const bool fixed_status )
2345 const int max_iter = this->param_.max_inner_iter_wells_;
2347 bool converged =
false;
2348 bool relax_convergence =
false;
2349 this->regularize_ =
false;
2350 const auto& summary_state = simulator.vanguard().summaryState();
2355 constexpr int min_its_after_switch = 4;
2356 int its_since_last_switch = min_its_after_switch;
2357 int switch_count= 0;
2359 const auto well_status_orig = this->wellStatus_;
2360 const auto operability_orig = this->operability_status_;
2361 auto well_status_cur = well_status_orig;
2362 int status_switch_count = 0;
2364 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2365 well_state.
well(this->index_of_well_).status == WellStatus::OPEN;
2367 const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) &&
2368 (!fixed_control || !fixed_status) && allow_open;
2369 bool changed =
false;
2370 bool final_check =
false;
2372 this->operability_status_.resetOperability();
2373 this->operability_status_.solvable =
true;
2375 its_since_last_switch++;
2376 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2377 const double wqTotal = this->primary_variables_.eval(WQTotal).value();
2378 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger, fixed_control, fixed_status);
2380 its_since_last_switch = 0;
2382 if (well_status_cur != this->wellStatus_) {
2383 well_status_cur = this->wellStatus_;
2384 status_switch_count++;
2387 if (!changed && final_check) {
2390 final_check =
false;
2394 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2396 if (it > this->param_.strict_inner_iter_wells_) {
2397 relax_convergence =
true;
2398 this->regularize_ =
true;
2401 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2403 converged = report.converged();
2407 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2409 its_since_last_switch = min_its_after_switch;
2416 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2417 initPrimaryVariablesEvaluation();
2419 }
while (it < max_iter);
2422 if (allow_switching){
2424 const bool is_stopped = this->wellIsStopped();
2425 if (this->wellHasTHPConstraints(summary_state)){
2426 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2427 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2429 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2433 this->wellStatus_ = well_status_orig;
2434 this->operability_status_ = operability_orig;
2435 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2436 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2437 deferred_logger.
debug(message);
2443 template<
typename TypeTag>
2450 std::vector<double> well_q_s(this->num_components_, 0.);
2451 const EvalWell&
bhp = this->primary_variables_.eval(Bhp);
2452 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2453 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2454 const int cell_idx = this->well_cells_[perf];
2455 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2456 std::vector<Scalar> mob(this->num_components_, 0.);
2457 getMobility(simulator, perf, mob, deferred_logger);
2458 std::vector<Scalar> cq_s(this->num_components_, 0.);
2459 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
2460 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2461 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2463 computePerfRate(intQuants, mob,
bhp.
value(), Tw, perf, allow_cf,
2464 cq_s, perf_rates, deferred_logger);
2465 for (
int comp = 0; comp < this->num_components_; ++comp) {
2466 well_q_s[comp] += cq_s[comp];
2469 const auto& comm = this->parallel_well_info_.communication();
2470 if (comm.size() > 1)
2472 comm.sum(well_q_s.data(), well_q_s.size());
2479 template <
typename TypeTag>
2484 const int num_pri_vars = this->primary_variables_.numWellEq();
2485 std::vector<double> retval(num_pri_vars);
2486 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2487 retval[ii] = this->primary_variables_.value(ii);
2496 template <
typename TypeTag>
2501 const int num_pri_vars = this->primary_variables_.numWellEq();
2502 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2503 this->primary_variables_.setValue(ii, it[ii]);
2505 return num_pri_vars;
2509 template <
typename TypeTag>
2513 const std::vector<EvalWell>& cq_s,
2514 const IntensiveQuantities& intQuants,
2517 auto fs = intQuants.fluidState();
2519 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2520 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2525 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2526 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2527 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2528 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2529 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2532 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2533 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2538 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2540 deferred_logger.
debug(
2541 fmt::format(
"Problematic d value {} obtained for well {}"
2542 " during calculateSinglePerf with rs {}"
2543 ", rv {}. Continue as if no dissolution (rs = 0) and"
2544 " vaporization (rv = 0) for this connection.",
2545 d, this->name(), fs.Rs(), fs.Rv()));
2546 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2548 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2549 cq_r_thermal = (cq_s[gasCompIdx] -
2550 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2551 (d * this->extendEval(fs.invB(phaseIdx)) );
2552 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2554 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2556 (d * this->extendEval(fs.invB(phaseIdx)) );
2562 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
2564 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2565 fs.setTemperature(this->well_ecl_.temperature());
2566 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2567 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2568 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2569 paramCache.setRegionIndex(pvtRegionIdx);
2570 paramCache.setMaxOilSat(maxOilSaturation);
2571 paramCache.updatePhase(fs, phaseIdx);
2573 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2574 fs.setDensity(phaseIdx, rho);
2575 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2576 fs.setEnthalpy(phaseIdx, h);
2577 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2578 result += getValue(cq_r_thermal);
2581 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2582 result += Base::restrictEval(cq_r_thermal);
2590 template <
typename TypeTag>
2591 template<
class Value>
2593 StandardWell<TypeTag>::
2594 gasOilPerfRateInj(
const std::vector<Value>& cq_s,
2595 PerforationRates& perf_rates,
2598 const Value& pressure,
2600 DeferredLogger& deferred_logger)
const
2602 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2603 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2612 const double d = 1.0 - getValue(rv) * getValue(rs);
2615 deferred_logger.debug(dValueError(d, this->name(),
2616 "gasOilPerfRateInj",
2621 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2624 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
2627 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2632 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2638 template <
typename TypeTag>
2639 template<
class Value>
2641 StandardWell<TypeTag>::
2642 gasOilPerfRateProd(std::vector<Value>& cq_s,
2643 PerforationRates& perf_rates,
2646 const Value& rvw)
const
2648 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2649 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2650 const Value cq_sOil = cq_s[oilCompIdx];
2651 const Value cq_sGas = cq_s[gasCompIdx];
2652 const Value dis_gas = rs * cq_sOil;
2653 const Value vap_oil = rv * cq_sGas;
2655 cq_s[gasCompIdx] += dis_gas;
2656 cq_s[oilCompIdx] += vap_oil;
2659 if (this->isProducer()) {
2660 perf_rates.dis_gas = getValue(dis_gas);
2661 perf_rates.vap_oil = getValue(vap_oil);
2664 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2665 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2666 const Value vap_wat = rvw * cq_sGas;
2667 cq_s[waterCompIdx] += vap_wat;
2668 if (this->isProducer())
2669 perf_rates.vap_wat = getValue(vap_wat);
2674 template <
typename TypeTag>
2675 template<
class Value>
2677 StandardWell<TypeTag>::
2678 gasWaterPerfRateProd(std::vector<Value>& cq_s,
2679 PerforationRates& perf_rates,
2681 const Value& rsw)
const
2683 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2684 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2685 const Value cq_sWat = cq_s[waterCompIdx];
2686 const Value cq_sGas = cq_s[gasCompIdx];
2687 const Value vap_wat = rvw * cq_sGas;
2688 const Value dis_gas_wat = rsw * cq_sWat;
2689 cq_s[waterCompIdx] += vap_wat;
2690 cq_s[gasCompIdx] += dis_gas_wat;
2691 if (this->isProducer()) {
2692 perf_rates.vap_wat = getValue(vap_wat);
2693 perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
2698 template <
typename TypeTag>
2699 template<
class Value>
2701 StandardWell<TypeTag>::
2702 gasWaterPerfRateInj(
const std::vector<Value>& cq_s,
2703 PerforationRates& perf_rates,
2706 const Value& pressure,
2707 DeferredLogger& deferred_logger)
const
2710 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2711 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2713 const double dw = 1.0 - getValue(rvw) * getValue(rsw);
2716 deferred_logger.debug(dValueError(dw, this->name(),
2717 "gasWaterPerfRateInj",
2718 rsw, rvw, pressure));
2722 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
2725 perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
2730 template <
typename TypeTag>
2731 template<
class Value>
2733 StandardWell<TypeTag>::
2734 disOilVapWatVolumeRatio(Value& volumeRatio,
2737 const Value& pressure,
2738 const std::vector<Value>& cmix_s,
2739 const std::vector<Value>& b_perfcells_dense,
2740 DeferredLogger& deferred_logger)
const
2742 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2743 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2745 const Value d = 1.0 - rvw * rsw;
2748 deferred_logger.debug(dValueError(d, this->name(),
2749 "disOilVapWatVolumeRatio",
2750 rsw, rvw, pressure));
2752 const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
2753 : cmix_s[waterCompIdx];
2754 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
2756 const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
2757 : cmix_s[waterCompIdx];
2758 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
2762 template <
typename TypeTag>
2763 template<
class Value>
2765 StandardWell<TypeTag>::
2766 gasOilVolumeRatio(Value& volumeRatio,
2769 const Value& pressure,
2770 const std::vector<Value>& cmix_s,
2771 const std::vector<Value>& b_perfcells_dense,
2772 DeferredLogger& deferred_logger)
const
2774 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2775 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2777 const Value d = 1.0 - rv * rs;
2780 deferred_logger.debug(dValueError(d, this->name(),
2781 "gasOilVolumeRatio",
2784 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
2785 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
2787 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
2788 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
#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:35
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:184
Class handling assemble of the equation system for StandardWell.
Definition: StandardWellAssemble.hpp:43
PrimaryVariables primary_variables_
Primary variables for well.
Definition: StandardWellEval.hpp:101
Definition: StandardWell.hpp:59
void computeWellConnectionPressures(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1340
typename StdWellEval::EvalWell EvalWell
Definition: StandardWell.hpp:120
typename StdWellEval::BVectorWell BVectorWell
Definition: StandardWell.hpp:121
void checkOperabilityUnderBHPLimit(const WellState< Scalar > &well_state, const Simulator &simulator, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:925
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:2288
std::vector< double > getPrimaryVars() const override
Definition: StandardWell_impl.hpp:2482
std::optional< double > computeBhpAtThpLimitProdWithAlq(const Simulator &simulator, const SummaryState &summary_state, const double alq_value, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2178
double connectionDensity(const int globalConnIdx, const int openConnIdx) const override
Definition: StandardWell_impl.hpp:1760
void updatePrimaryVariables(const SummaryState &summary_state, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1775
void addWellContributions(SparseMatrixAdapter &mat) const override
Definition: StandardWell_impl.hpp:1873
void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator &wellPICalc, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1222
EvalWell pskinwater(const double throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1901
void initPrimaryVariablesEvaluation() override
Definition: StandardWell_impl.hpp:113
void checkConvergenceExtraEqs(const std::vector< double > &res, ConvergenceReport &report) const
Definition: StandardWell_impl.hpp:2102
void updateIPRImplicit(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:849
std::vector< double > computeWellPotentialWithTHP(const Simulator &simulator, DeferredLogger &deferred_logger, const WellState< Scalar > &well_state) const
Definition: StandardWell_impl.hpp:1561
void computePropertiesForWellConnectionPressures(const Simulator &simulator, const WellState< Scalar > &well_state, WellConnectionProps &props) const
Definition: StandardWell_impl.hpp:1132
EvalWell pskin(const double throuhgput, const EvalWell &water_velocity, const EvalWell &poly_inj_conc, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1933
void updateWaterMobilityWithPolymer(const Simulator &simulator, const int perf, std::vector< EvalWell > &mob_water, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1809
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:2335
void updateWaterThroughput(const double dt, WellState< Scalar > &well_state) const override
Definition: StandardWell_impl.hpp:2007
int setPrimaryVars(std::vector< double >::const_iterator it) override
Definition: StandardWell_impl.hpp:2499
std::optional< double > computeBhpAtThpLimitInj(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2256
void updatePrimaryVariablesNewton(const BVectorWell &dwells, const bool stop_or_zero_rate_target, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:716
void calculateExplicitQuantities(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1382
void computeWellRatesWithThpAlqProd(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< double > &potentials, double alq) const
Definition: StandardWell_impl.hpp:1686
std::vector< double > computeCurrentWellRates(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2446
StandardWell(const Well &well, const ParallelWellInfo &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 > &perf_data)
Definition: StandardWell_impl.hpp:72
double computeWellRatesAndBhpWithThpAlqProd(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< double > &potentials, double alq) const
Definition: StandardWell_impl.hpp:1658
typename StdWellEval::StdWellConnections::Properties WellConnectionProps
Definition: StandardWell.hpp:273
double getRefDensity() const override
Definition: StandardWell_impl.hpp:1798
void updateWellState(const SummaryState &summary_state, const BVectorWell &dwells, WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:695
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:639
virtual void init(const PhaseUsage *phase_usage_arg, const std::vector< double > &depth_arg, const double gravity_arg, const int num_cells, const std::vector< Scalar > &B_avg, const bool changed_to_open_this_step) override
Definition: StandardWell_impl.hpp:96
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 &perf_rates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:126
void updateIPR(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:759
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:338
void computeWellConnectionDensitesPressures(const Simulator &simulator, const WellState< Scalar > &well_state, const WellConnectionProps &props, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1298
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1398
typename StdWellEval::Eval Eval
Definition: StandardWell.hpp:119
void computeWellPotentials(const Simulator &simulator, const WellState< Scalar > &well_state, std::vector< double > &well_potentials, DeferredLogger &deferred_logger) override
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1703
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:2121
bool openCrossFlowAvoidSingularity(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1121
std::optional< double > computeBhpAtThpLimitProd(const WellState< Scalar > &well_state, const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2164
bool allDrawDownWrongDirection(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1043
void recoverWellSolutionAndUpdateWellState(const SummaryState &summary_state, const BVector &x, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1430
void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target, WellState< Scalar > &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:739
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:2054
virtual ConvergenceReport getWellConvergence(const SummaryState &summary_state, const WellState< Scalar > &well_state, const std::vector< double > &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:1178
void computeWellRatesWithBhp(const Simulator &simulator, const double &bhp, std::vector< double > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1450
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:364
void checkOperabilityUnderTHPLimit(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:994
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:1881
bool canProduceInjectWithCurrentBhp(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1085
void handleInjectivityRate(const Simulator &simulator, const int perf, std::vector< EvalWell > &cq_s) const
Definition: StandardWell_impl.hpp:2031
EvalWell wpolymermw(const double throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1978
bool computeWellPotentialsImplicit(const Simulator &simulator, std::vector< double > &well_potentials, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1595
void solveEqAndUpdateWellState(const SummaryState &summary_state, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1360
void computeWellRatesWithBhpIterations(const Simulator &simulator, const double &bhp, std::vector< double > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1496
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:489
EvalWell getQs(const int compIdx) const
Returns scaled rate for a component.
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:42
std::optional< double > computeBhpAtThpLimitProd(const std::function< std::vector< double >(const double)> &frates, const SummaryState &summary_state, const double maxPerfPress, const double rho, const double alq_value, const double thp_limit, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for a producer.
double mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
double calculateThpFromBhp(const std::vector< double > &rates, const double bhp, const double rho, const std::optional< double > &alq, const double thp_limit, DeferredLogger &deferred_logger) const
Calculates THP from BHP.
std::optional< double > computeBhpAtThpLimitInj(const std::function< std::vector< double >(const double)> &frates, const SummaryState &summary_state, const double rho, const double 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:37
std::pair< bool, bool > computeWellPotentials(std::vector< double > &well_potentials, const WellState< double > &well_state)
void prepareForPotentialCalculations(const SummaryState &summary_state, WellState< double > &well_state, Well::InjectionControls &inj_controls, Well::ProductionControls &prod_controls) const
const int num_components_
Definition: WellInterfaceGeneric.hpp:311
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:302
Definition: WellInterfaceIndices.hpp:35
Definition: WellInterface.hpp:75
Dune::BCRSMatrix< Opm::MatrixBlock< double, 1, 1 > > PressureMatrix
Definition: WellInterface.hpp:102
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:82
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1691
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:105
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:85
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:96
Dune::BlockVector< VectorBlockType > BVector
Definition: WellInterface.hpp:101
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:83
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:84
bool solveWellWithTHPConstraint(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:486
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: WellInterface.hpp:87
Definition: WellProdIndexCalculator.hpp:36
double connectionProdIndStandard(const std::size_t connIdx, const double connMobility) const
Definition: WellState.hpp:62
const SingleWellState< Scalar > & well(std::size_t well_index) const
Definition: WellState.hpp:300
std::vector< Scalar > & wellRates(std::size_t well_index)
One rate per well and phase.
Definition: WellState.hpp:265
@ NONE
Definition: DeferredLogger.hpp:46
VFPEvaluation bhp(const VFPProdTable &table, const double aqua, const double liquid, const double vapour, const double thp, const double alq, const double explicit_wfr, const double explicit_gfr, const bool use_vfpexplicit)
Definition: BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:484
Definition: PerforationData.hpp:40
double dis_gas
Definition: PerforationData.hpp:41
double vap_oil
Definition: PerforationData.hpp:43
double vap_wat
Definition: PerforationData.hpp:44
double dis_gas_in_water
Definition: PerforationData.hpp:42
Definition: BlackoilPhases.hpp:46
double value
Definition: VFPHelpers.hpp:105