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,
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<Scalar>& depth_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);
289 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
290 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
291 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
294 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
295 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
296 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
299 if constexpr (Indices::enableSolvent) {
300 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
303 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
304 gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
305 cmix_s, b_perfcells_dense, deferred_logger);
307 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
308 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
309 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
311 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
312 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
313 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
319 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
320 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
321 Value cqt_is = cqt_i / volumeRatio;
322 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
326 if (this->isProducer()) {
327 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
328 gasOilPerfRateInj(cq_s, perf_rates,
329 rv, rs, pressure, rvw, deferred_logger);
331 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
333 gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
334 pressure, deferred_logger);
341 template<
typename TypeTag>
346 const Well::InjectionControls& inj_controls,
347 const Well::ProductionControls& prod_controls,
354 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
357 this->linSys_.clear();
359 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
360 prod_controls, well_state,
361 group_state, deferred_logger);
367 template<
typename TypeTag>
372 const Well::InjectionControls& inj_controls,
373 const Well::ProductionControls& prod_controls,
379 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
380 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
382 auto& ws = well_state.
well(this->index_of_well_);
383 ws.phase_mixing_rates.fill(0.0);
386 const int np = this->number_of_phases_;
388 std::vector<RateVector> connectionRates = this->connectionRates_;
390 auto& perf_data = ws.perf_data;
391 auto& perf_rates = perf_data.phase_rates;
392 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
394 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
395 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
396 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
397 calculateSinglePerf(simulator, perf, well_state, connectionRates,
398 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
401 if constexpr (has_polymer && Base::has_polymermw) {
402 if (this->isInjector()) {
403 handleInjectivityEquations(simulator, well_state, perf,
404 water_flux_s, deferred_logger);
407 const int cell_idx = this->well_cells_[perf];
408 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
410 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
412 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
415 assemblePerforationEq(cq_s_effective,
418 this->primary_variables_.numWellEq(),
422 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
423 auto& perf_rate_solvent = perf_data.solvent_rates;
424 perf_rate_solvent[perf] = cq_s[componentIdx].value();
426 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
430 if constexpr (has_zFraction) {
432 assembleZFracEq(cq_s_zfrac_effective,
434 this->primary_variables_.numWellEq(),
439 this->connectionRates_ = connectionRates;
444 const auto& comm = this->parallel_well_info_.communication();
445 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
449 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
452 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
455 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
456 if (FluidSystem::numActivePhases() > 1) {
458 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
459 this->F0_[componentIdx]) * volume / dt;
461 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
463 assembleSourceEq(resWell_loc,
465 this->primary_variables_.numWellEq(),
469 const auto& summaryState = simulator.vanguard().summaryState();
470 const Schedule& schedule = simulator.vanguard().schedule();
471 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
473 assembleControlEq(well_state, group_state,
474 schedule, summaryState,
475 inj_controls, prod_controls,
476 this->primary_variables_,
477 this->connections_.rho(),
479 stopped_or_zero_target,
485 this->linSys_.invert();
487 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
494 template<
typename TypeTag>
500 std::vector<RateVector>& connectionRates,
501 std::vector<EvalWell>& cq_s,
506 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
507 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
508 const int cell_idx = this->well_cells_[perf];
509 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
510 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
511 getMobility(simulator, perf, mob, deferred_logger);
514 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
515 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
516 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
517 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
518 cq_s, perf_rates, deferred_logger);
520 auto& ws = well_state.
well(this->index_of_well_);
521 auto& perf_data = ws.perf_data;
522 if constexpr (has_polymer && Base::has_polymermw) {
523 if (this->isInjector()) {
526 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
527 water_flux_s = cq_s[water_comp_idx];
530 handleInjectivityRate(simulator, perf, cq_s);
535 if (this->isProducer()) {
536 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.
dis_gas;
537 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.
dis_gas_in_water;
538 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.
vap_oil;
539 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.
vap_wat;
540 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.
dis_gas;
541 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.
dis_gas_in_water;
542 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.
vap_oil;
543 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.
vap_wat;
546 if constexpr (has_energy) {
547 connectionRates[perf][Indices::contiEnergyEqIdx] =
548 connectionRateEnergy(simulator.problem().maxOilSaturation(cell_idx),
549 cq_s, intQuants, deferred_logger);
552 if constexpr (has_polymer) {
553 std::variant<Scalar,EvalWell> polymerConcentration;
554 if (this->isInjector()) {
555 polymerConcentration = this->wpolymer();
557 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
558 intQuants.polymerViscosityCorrection());
561 [[maybe_unused]]
EvalWell cq_s_poly;
562 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
564 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
565 cq_s, polymerConcentration);
567 if constexpr (Base::has_polymermw) {
568 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
569 perf, connectionRates, deferred_logger);
573 if constexpr (has_foam) {
574 std::variant<Scalar,EvalWell> foamConcentration;
575 if (this->isInjector()) {
576 foamConcentration = this->wfoam();
578 foamConcentration = this->extendEval(intQuants.foamConcentration());
580 connectionRates[perf][Indices::contiFoamEqIdx] =
581 this->connections_.connectionRateFoam(cq_s, foamConcentration,
582 FoamModule::transportPhase(),
586 if constexpr (has_zFraction) {
587 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
588 if (this->isInjector()) {
589 solventConcentration = this->wsolvent();
591 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
592 this->extendEval(intQuants.yVolume())};
594 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
595 cq_s_zfrac_effective) =
596 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
598 solventConcentration);
601 if constexpr (has_brine) {
602 std::variant<Scalar,EvalWell> saltConcentration;
603 if (this->isInjector()) {
604 saltConcentration = this->wsalt();
606 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
609 connectionRates[perf][Indices::contiBrineEqIdx] =
610 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
615 if constexpr (has_micp) {
616 std::variant<Scalar,EvalWell> microbialConcentration;
617 std::variant<Scalar,EvalWell> oxygenConcentration;
618 std::variant<Scalar,EvalWell> ureaConcentration;
619 if (this->isInjector()) {
620 microbialConcentration = this->wmicrobes();
621 oxygenConcentration = this->woxygen();
622 ureaConcentration = this->wurea();
624 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
625 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
626 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
628 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
629 connectionRates[perf][Indices::contiOxygenEqIdx],
630 connectionRates[perf][Indices::contiUreaEqIdx]) =
631 this->connections_.connectionRatesMICP(cq_s,
632 microbialConcentration,
638 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
643 template<
typename TypeTag>
644 template<
class Value>
649 std::vector<Value>& mob,
652 auto obtain = [
this](
const Eval& value)
654 if constexpr (std::is_same_v<Value, Scalar>) {
655 static_cast<void>(
this);
656 return getValue(value);
658 return this->extendEval(value);
662 obtain, deferred_logger);
665 if constexpr (has_polymer) {
666 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
667 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
672 if constexpr (!Base::has_polymermw) {
673 if constexpr (std::is_same_v<Value, Scalar>) {
674 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
675 for (std::size_t i = 0; i < mob.size(); ++i) {
676 mob_eval[i].setValue(mob[i]);
678 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
679 for (std::size_t i = 0; i < mob.size(); ++i) {
680 mob[i] = getValue(mob_eval[i]);
683 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
690 const Scalar bhp = this->primary_variables_.value(Bhp);
691 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
692 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press);
693 for (std::size_t i = 0; i < mob.size(); ++i) {
694 mob[i] *= multiplier;
700 template<
typename TypeTag>
708 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
710 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
711 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
713 const auto& summary_state = simulator.vanguard().summaryState();
714 updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
715 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.
well(this->index_of_well_));
722 template<
typename TypeTag>
726 const bool stop_or_zero_rate_target,
729 const Scalar dFLimit = this->param_.dwell_fraction_max_;
730 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
731 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
734 if constexpr (Base::has_polymermw) {
735 this->primary_variables_.updateNewtonPolyMW(dwells);
738 this->primary_variables_.checkFinite(deferred_logger);
745 template<
typename TypeTag>
750 const SummaryState& summary_state,
753 this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
756 if constexpr (Base::has_polymermw) {
757 this->primary_variables_.copyToWellStatePolyMW(well_state);
765 template<
typename TypeTag>
773 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
774 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
776 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
777 std::vector<Scalar> mob(this->num_components_, 0.0);
778 getMobility(simulator, perf, mob, deferred_logger);
780 const int cell_idx = this->well_cells_[perf];
781 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
782 const auto& fs = int_quantities.fluidState();
784 Scalar p_r = this->getPerfCellPressure(fs).value();
787 std::vector<Scalar> b_perf(this->num_components_);
788 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
789 if (!FluidSystem::phaseIsActive(phase)) {
792 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
793 b_perf[comp_idx] = fs.invB(phase).value();
795 if constexpr (has_solvent) {
796 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
800 const Scalar h_perf = this->connections_.pressure_diff(perf);
801 const Scalar pressure_diff = p_r - h_perf;
806 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
807 deferred_logger.
debug(
"CROSSFLOW_IPR",
808 "cross flow found when updateIPR for well " + name()
809 +
" . The connection is ignored in IPR calculations");
815 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
816 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
817 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
821 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
822 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
823 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
824 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
825 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
826 ipr_b_perf[comp_idx] += tw_mob;
830 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
831 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
832 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
833 const Scalar rs = (fs.Rs()).value();
834 const Scalar rv = (fs.Rv()).value();
836 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
837 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
839 ipr_a_perf[gas_comp_idx] += dis_gas_a;
840 ipr_a_perf[oil_comp_idx] += vap_oil_a;
842 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
843 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
845 ipr_b_perf[gas_comp_idx] += dis_gas_b;
846 ipr_b_perf[oil_comp_idx] += vap_oil_b;
849 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
850 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
851 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
854 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
855 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
858 template<
typename TypeTag>
872 auto rates = well_state.
well(this->index_of_well_).surface_rates;
874 for (std::size_t p = 0; p < rates.size(); ++p) {
875 zero_rates &= rates[p] == 0.0;
877 auto& ws = well_state.
well(this->index_of_well_);
879 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
880 deferred_logger.
debug(msg);
892 const auto& group_state = simulator.problem().wellModel().groupState();
894 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
895 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
897 auto inj_controls = Well::InjectionControls(0);
898 auto prod_controls = Well::ProductionControls(0);
899 prod_controls.addControl(Well::ProducerCMode::BHP);
900 prod_controls.bhp_limit = well_state.
well(this->index_of_well_).bhp;
903 const auto cmode = ws.production_cmode;
904 ws.production_cmode = Well::ProducerCMode::BHP;
905 const double dt = simulator.timeStepSize();
906 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
908 const size_t nEq = this->primary_variables_.numWellEq();
912 for (
size_t i=0; i < nEq; ++i){
918 x_well[0].resize(nEq);
919 this->linSys_.solve(rhs, x_well);
921 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
922 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
923 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
924 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
926 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
928 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
931 ws.production_cmode = cmode;
934 template<
typename TypeTag>
941 const auto& summaryState = simulator.vanguard().summaryState();
945 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
946 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
949 Scalar total_ipr_mass_rate = 0.0;
950 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
952 if (!FluidSystem::phaseIsActive(phaseIdx)) {
956 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
957 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
959 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
960 total_ipr_mass_rate += ipr_rate * rho;
962 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
963 this->operability_status_.operable_under_only_bhp_limit =
false;
967 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
971 std::vector<Scalar> well_rates_bhp_limit;
972 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
974 this->adaptRatesForVFP(well_rates_bhp_limit);
975 const Scalar thp_limit = this->getTHPConstraint(summaryState);
978 this->connections_.rho(),
979 this->getALQ(well_state),
982 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
983 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
994 this->operability_status_.operable_under_only_bhp_limit =
true;
995 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1003 template<
typename TypeTag>
1010 const auto& summaryState = simulator.vanguard().summaryState();
1011 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1012 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1015 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1018 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1019 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1021 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1022 if (this->isProducer() && *obtain_bhp < thp_limit) {
1023 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1024 +
" bars is SMALLER than thp limit "
1026 +
" bars as a producer for well " + name();
1027 deferred_logger.
debug(msg);
1029 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1030 const std::string msg =
" obtained bhp " +
std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1031 +
" bars is LARGER than thp limit "
1033 +
" bars as a injector for well " + name();
1034 deferred_logger.
debug(msg);
1037 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1038 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1039 if (!this->wellIsStopped()) {
1040 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1041 deferred_logger.
debug(
" could not find bhp value at thp limit "
1043 +
" bar for well " + name() +
", the well might need to be closed ");
1052 template<
typename TypeTag>
1057 bool all_drawdown_wrong_direction =
true;
1059 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1060 const int cell_idx = this->well_cells_[perf];
1061 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1062 const auto& fs = intQuants.fluidState();
1064 const Scalar pressure = this->getPerfCellPressure(fs).value();
1065 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1068 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1069 const Scalar drawdown = pressure - well_pressure;
1074 if ( (drawdown < 0. && this->isInjector()) ||
1075 (drawdown > 0. && this->isProducer()) ) {
1076 all_drawdown_wrong_direction =
false;
1081 const auto& comm = this->parallel_well_info_.communication();
1082 if (comm.size() > 1)
1084 all_drawdown_wrong_direction =
1085 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1088 return all_drawdown_wrong_direction;
1094 template<
typename TypeTag>
1101 const Scalar bhp = well_state.
well(this->index_of_well_).bhp;
1102 std::vector<Scalar> well_rates;
1103 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1105 const Scalar sign = (this->isProducer()) ? -1. : 1.;
1106 const Scalar threshold = sign * std::numeric_limits<Scalar>::min();
1108 bool can_produce_inject =
false;
1109 for (
const auto value : well_rates) {
1110 if (this->isProducer() && value < threshold) {
1111 can_produce_inject =
true;
1113 }
else if (this->isInjector() && value > threshold) {
1114 can_produce_inject =
true;
1119 if (!can_produce_inject) {
1120 deferred_logger.
debug(
" well " + name() +
" CANNOT produce or inejct ");
1123 return can_produce_inject;
1130 template<
typename TypeTag>
1135 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1141 template<
typename TypeTag>
1147 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1149 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1151 return model.intensiveQuantities(cell_idx, 0)
1152 .fluidState().temperature(phase_idx).value();
1156 [&model = simulator.model()](
int cell_idx)
1158 return model.intensiveQuantities(cell_idx, 0)
1159 .fluidState().saltConcentration().value();
1163 [&model = simulator.model()](
int cell_idx)
1165 return model.intensiveQuantities(cell_idx, 0)
1166 .fluidState().pvtRegionIndex();
1170 if constexpr (Indices::enableSolvent) {
1171 prop_func.solventInverseFormationVolumeFactor =
1172 [&model = simulator.model()](
int cell_idx)
1174 return model.intensiveQuantities(cell_idx, 0)
1175 .solventInverseFormationVolumeFactor().value();
1178 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1180 return model.intensiveQuantities(cell_idx, 0)
1181 .solventRefDensity();
1185 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1192 template<
typename TypeTag>
1197 const std::vector<Scalar>& B_avg,
1199 const bool relax_tolerance)
const
1203 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1205 Scalar tol_wells = this->param_.tolerance_wells_;
1207 constexpr Scalar stopped_factor = 1.e-4;
1209 constexpr Scalar dynamic_thp_factor = 1.e-1;
1210 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1211 tol_wells = tol_wells*stopped_factor;
1212 }
else if (this->getDynamicThpLimit()) {
1213 tol_wells = tol_wells*dynamic_thp_factor;
1216 std::vector<Scalar> res;
1219 this->param_.max_residual_allowed_,
1221 this->param_.relaxed_tolerance_flow_well_,
1223 this->wellIsStopped(),
1227 checkConvergenceExtraEqs(res, report);
1236 template<
typename TypeTag>
1244 auto fluidState = [&simulator,
this](
const int perf)
1246 const auto cell_idx = this->well_cells_[perf];
1247 return simulator.model()
1248 .intensiveQuantities(cell_idx, 0).fluidState();
1251 const int np = this->number_of_phases_;
1252 auto setToZero = [np](
Scalar* x) ->
void
1254 std::fill_n(x, np, 0.0);
1257 auto addVector = [np](
const Scalar* src,
Scalar* dest) ->
void
1259 std::transform(src, src + np, dest, dest, std::plus<>{});
1262 auto& ws = well_state.
well(this->index_of_well_);
1263 auto& perf_data = ws.perf_data;
1264 auto* wellPI = ws.productivity_index.data();
1265 auto* connPI = perf_data.prod_index.data();
1269 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1270 auto subsetPerfID = 0;
1272 for (
const auto& perf : *this->perf_data_) {
1273 auto allPerfID = perf.ecl_index;
1275 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) ->
Scalar
1280 std::vector<Scalar> mob(this->num_components_, 0.0);
1281 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1283 const auto& fs = fluidState(subsetPerfID);
1286 if (this->isInjector()) {
1287 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1288 mob, connPI, deferred_logger);
1291 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1294 addVector(connPI, wellPI);
1301 const auto& comm = this->parallel_well_info_.communication();
1302 if (comm.size() > 1) {
1303 comm.sum(wellPI, np);
1306 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1307 "Internal logic error in processing connections for PI/II");
1312 template<
typename TypeTag>
1322 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1327 [&model = simulator.model()](
const int cell,
1328 const std::vector<int>& phases,
1329 std::vector<Scalar>& mob)
1331 const auto& iq = model.intensiveQuantities(cell, 0);
1333 std::transform(phases.begin(), phases.end(), mob.begin(),
1334 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1339 [&model = simulator.model()](
const int cell,
1340 const std::vector<int>& phases,
1341 std::vector<Scalar>& rho)
1343 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1345 std::transform(phases.begin(), phases.end(), rho.begin(),
1346 [&fs](
const int phase) { return fs.density(phase).value(); });
1350 const auto stopped_or_zero_rate_target = this->
1351 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1354 .computeProperties(stopped_or_zero_rate_target, well_state,
1355 prop_func, props, deferred_logger);
1362 template<
typename TypeTag>
1369 const auto props = computePropertiesForWellConnectionPressures
1370 (simulator, well_state);
1372 computeWellConnectionDensitesPressures(simulator, well_state,
1373 props, deferred_logger);
1380 template<
typename TypeTag>
1387 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1392 dx_well[0].resize(this->primary_variables_.numWellEq());
1393 this->linSys_.solve( dx_well);
1395 updateWellState(simulator, dx_well, well_state, deferred_logger);
1402 template<
typename TypeTag>
1409 updatePrimaryVariables(simulator, well_state, deferred_logger);
1410 initPrimaryVariablesEvaluation();
1411 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1412 this->computeAccumWell();
1417 template<
typename TypeTag>
1422 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1424 if (this->param_.matrix_add_well_contributions_)
1430 this->linSys_.apply(x, Ax);
1436 template<
typename TypeTag>
1441 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1443 this->linSys_.apply(r);
1449 template<
typename TypeTag>
1457 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1460 xw[0].resize(this->primary_variables_.numWellEq());
1462 this->linSys_.recoverSolutionWell(x, xw);
1463 updateWellState(simulator, xw, well_state, deferred_logger);
1469 template<
typename TypeTag>
1474 std::vector<Scalar>& well_flux,
1478 const int np = this->number_of_phases_;
1479 well_flux.resize(np, 0.0);
1481 const bool allow_cf = this->getAllowCrossFlow();
1483 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1484 const int cell_idx = this->well_cells_[perf];
1485 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1487 std::vector<Scalar> mob(this->num_components_, 0.);
1488 getMobility(simulator, perf, mob, deferred_logger);
1489 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1490 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1491 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1493 std::vector<Scalar> cq_s(this->num_components_, 0.);
1495 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1496 cq_s, perf_rates, deferred_logger);
1498 for(
int p = 0; p < np; ++p) {
1499 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1503 if constexpr (has_solvent) {
1505 assert(pu.phase_used[Gas]);
1506 const int gas_pos = pu.phase_pos[Gas];
1507 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1510 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1515 template<
typename TypeTag>
1520 std::vector<Scalar>& well_flux,
1530 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1531 const auto& group_state = simulator.problem().wellModel().groupState();
1534 const auto& summary_state = simulator.vanguard().summaryState();
1535 auto inj_controls = well_copy.
well_ecl_.isInjector()
1536 ? well_copy.
well_ecl_.injectionControls(summary_state)
1537 : Well::InjectionControls(0);
1538 auto prod_controls = well_copy.
well_ecl_.isProducer()
1539 ? well_copy.
well_ecl_.productionControls(summary_state) :
1540 Well::ProductionControls(0);
1543 auto& ws = well_state_copy.
well(this->index_of_well_);
1545 inj_controls.bhp_limit = bhp;
1546 ws.injection_cmode = Well::InjectorCMode::BHP;
1548 prod_controls.bhp_limit = bhp;
1549 ws.production_cmode = Well::ProducerCMode::BHP;
1554 const int np = this->number_of_phases_;
1555 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1556 for (
int phase = 0; phase < np; ++phase){
1557 well_state_copy.
wellRates(this->index_of_well_)[phase]
1558 = sign * ws.well_potentials[phase];
1564 const double dt = simulator.timeStepSize();
1565 const bool converged = well_copy.
iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1567 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1568 " potentials are computed based on unconverged solution";
1569 deferred_logger.
debug(msg);
1580 template<
typename TypeTag>
1581 std::vector<typename StandardWell<TypeTag>::Scalar>
1587 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1588 const auto& summary_state = simulator.vanguard().summaryState();
1590 const auto& well = this->well_ecl_;
1591 if (well.isInjector()){
1592 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1593 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1594 if (bhp_at_thp_limit) {
1595 const Scalar bhp = std::min(*bhp_at_thp_limit,
1596 static_cast<Scalar>(controls.bhp_limit));
1597 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1599 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1600 "Failed in getting converged thp based potential calculation for well "
1601 + name() +
". Instead the bhp based value is used");
1602 const Scalar bhp = controls.bhp_limit;
1603 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1606 computeWellRatesWithThpAlqProd(
1607 simulator, summary_state,
1608 deferred_logger, potentials, this->getALQ(well_state)
1615 template<
typename TypeTag>
1619 std::vector<Scalar>& well_potentials,
1628 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1629 const auto& group_state = simulator.problem().wellModel().groupState();
1630 auto& ws = well_state_copy.
well(this->index_of_well_);
1633 const auto& summary_state = simulator.vanguard().summaryState();
1634 auto inj_controls = well_copy.
well_ecl_.isInjector()
1635 ? well_copy.
well_ecl_.injectionControls(summary_state)
1636 : Well::InjectionControls(0);
1637 auto prod_controls = well_copy.
well_ecl_.isProducer()
1638 ? well_copy.
well_ecl_.productionControls(summary_state) :
1639 Well::ProductionControls(0);
1645 const int num_perf = ws.perf_data.size();
1646 for (
int perf = 0; perf < num_perf; ++perf) {
1650 const int np = this->number_of_phases_;
1651 bool trivial =
true;
1652 for (
int phase = 0; phase < np; ++phase){
1653 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1657 for (
int phase = 0; phase < np; ++phase) {
1658 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1663 const double dt = simulator.timeStepSize();
1665 bool converged =
false;
1666 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1667 converged = well_copy.
solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1669 converged = well_copy.
iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1673 well_potentials.clear();
1674 well_potentials.resize(np, 0.0);
1675 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1676 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1678 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1682 if constexpr (has_solvent) {
1684 assert(pu.phase_used[Gas]);
1685 const int gas_pos = pu.phase_pos[Gas];
1687 well_potentials[gas_pos] += rate.value();
1693 template<
typename TypeTag>
1697 const SummaryState &summary_state,
1699 std::vector<Scalar>& potentials,
1703 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1704 simulator, summary_state, alq, deferred_logger);
1705 if (bhp_at_thp_limit) {
1706 const auto& controls = this->well_ecl_.productionControls(summary_state);
1707 bhp = std::max(*bhp_at_thp_limit,
1708 static_cast<Scalar>(controls.bhp_limit));
1709 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1712 deferred_logger.
warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1713 "Failed in getting converged thp based potential calculation for well "
1714 + name() +
". Instead the bhp based value is used");
1715 const auto& controls = this->well_ecl_.productionControls(summary_state);
1716 bhp = controls.bhp_limit;
1717 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1722 template<
typename TypeTag>
1726 const SummaryState& summary_state,
1728 std::vector<Scalar>& potentials,
1732 computeWellRatesAndBhpWithThpAlqProd(simulator,
1739 template<
typename TypeTag>
1744 std::vector<Scalar>& well_potentials,
1747 const auto [compute_potential, bhp_controlled_well] =
1750 if (!compute_potential) {
1754 bool converged_implicit =
false;
1755 if (this->param_.local_well_solver_control_switching_) {
1756 converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
1758 if (!converged_implicit) {
1760 const auto& summaryState = simulator.vanguard().summaryState();
1761 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1771 const auto& ws = well_state.
well(this->index_of_well_);
1772 if (this->isInjector())
1773 bhp = std::max(ws.bhp, bhp);
1775 bhp = std::min(ws.bhp, bhp);
1777 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1778 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1781 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1785 this->checkNegativeWellPotentials(well_potentials,
1786 this->param_.check_well_operability_,
1796 template<
typename TypeTag>
1800 const int openConnIdx)
const
1802 return (openConnIdx < 0)
1804 : this->connections_.rho(openConnIdx);
1811 template<
typename TypeTag>
1818 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1820 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1821 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1824 if constexpr (Base::has_polymermw) {
1825 this->primary_variables_.updatePolyMW(well_state);
1828 this->primary_variables_.checkFinite(deferred_logger);
1834 template<
typename TypeTag>
1839 return this->connections_.rho();
1845 template<
typename TypeTag>
1850 std::vector<EvalWell>& mob,
1853 const int cell_idx = this->well_cells_[perf];
1854 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1855 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1859 if (this->isInjector()) {
1861 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1862 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1863 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1866 if (PolymerModule::hasPlyshlog()) {
1869 if (this->isInjector() && this->wpolymer() == 0.) {
1874 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1875 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1877 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1879 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1880 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1881 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1882 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1883 perf_rates, deferred_logger);
1885 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1886 const auto& material_law_manager = simulator.problem().materialLawManager();
1887 const auto& scaled_drainage_info =
1888 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1889 const Scalar swcr = scaled_drainage_info.Swcr;
1890 const EvalWell poro = this->extendEval(int_quant.porosity());
1891 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1893 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1894 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1895 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1897 if (PolymerModule::hasShrate()) {
1900 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1902 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1903 int_quant.pvtRegionIndex(),
1906 mob[waterCompIdx] /= shear_factor;
1910 template<
typename TypeTag>
1914 this->linSys_.extract(jacobian);
1918 template <
typename TypeTag>
1922 const int pressureVarIndex,
1923 const bool use_well_weights,
1926 this->linSys_.extractCPRPressureMatrix(jacobian,
1937 template<
typename TypeTag>
1944 if constexpr (Base::has_polymermw) {
1945 const int water_table_id = this->polymerWaterTable_();
1946 if (water_table_id <= 0) {
1948 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1951 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1952 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1954 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1955 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1959 fmt::format(
"Polymermw is not activated, while injecting "
1960 "skin pressure is requested for well {}", name()),
1969 template<
typename TypeTag>
1977 if constexpr (Base::has_polymermw) {
1978 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1979 const EvalWell water_velocity_abs = abs(water_velocity);
1980 if (poly_inj_conc == 0.) {
1981 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1983 const int polymer_table_id = this->polymerTable_();
1984 if (polymer_table_id <= 0) {
1986 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1989 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1990 const Scalar reference_concentration = skprpolytable.refConcentration;
1991 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1993 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1994 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1995 if (poly_inj_conc == reference_concentration) {
1996 return sign * pskin_poly;
1999 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2000 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2001 return sign * pskin;
2004 fmt::format(
"Polymermw is not activated, while injecting "
2005 "skin pressure is requested for well {}", name()),
2014 template<
typename TypeTag>
2021 if constexpr (Base::has_polymermw) {
2022 const int table_id = this->polymerInjTable_();
2023 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2024 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2025 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2026 if (this->wpolymer() == 0.) {
2027 return molecular_weight;
2029 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2030 return molecular_weight;
2033 fmt::format(
"Polymermw is not activated, while injecting "
2034 "polymer molecular weight is requested for well {}", name()),
2043 template<
typename TypeTag>
2049 if constexpr (Base::has_polymermw) {
2050 if (!this->isInjector()) {
2054 auto& perf_water_throughput = well_state.
well(this->index_of_well_)
2055 .perf_data.water_throughput;
2057 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2058 const Scalar perf_water_vel =
2059 this->primary_variables_.value(Bhp + 1 + perf);
2063 if (perf_water_vel >
Scalar{0}) {
2064 perf_water_throughput[perf] += perf_water_vel * dt;
2074 template<
typename TypeTag>
2079 std::vector<EvalWell>& cq_s)
const
2081 const int cell_idx = this->well_cells_[perf];
2082 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2083 const auto& fs = int_quants.fluidState();
2084 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2085 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2086 const int wat_vel_index = Bhp + 1 + perf;
2087 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2091 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2097 template<
typename TypeTag>
2106 const int cell_idx = this->well_cells_[perf];
2107 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2108 const auto& fs = int_quants.fluidState();
2109 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2110 const EvalWell water_flux_r = water_flux_s / b_w;
2111 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2112 const EvalWell water_velocity = water_flux_r / area;
2113 const int wat_vel_index = Bhp + 1 + perf;
2116 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2118 const auto& ws = well_state.
well(this->index_of_well_);
2119 const auto& perf_data = ws.perf_data;
2120 const auto& perf_water_throughput = perf_data.water_throughput;
2121 const Scalar throughput = perf_water_throughput[perf];
2122 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2124 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2125 poly_conc.setValue(this->wpolymer());
2128 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2129 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2132 assembleInjectivityEq(eq_pskin,
2137 this->primary_variables_.numWellEq(),
2145 template<
typename TypeTag>
2154 if constexpr (Base::has_polymermw) {
2156 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2164 template<
typename TypeTag>
2171 std::vector<RateVector>& connectionRates,
2176 if (this->isInjector()) {
2177 const int wat_vel_index = Bhp + 1 + perf;
2178 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2179 if (water_velocity > 0.) {
2180 const auto& ws = well_state.
well(this->index_of_well_);
2181 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2182 const Scalar throughput = perf_water_throughput[perf];
2183 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2184 cq_s_polymw *= molecular_weight;
2190 }
else if (this->isProducer()) {
2191 if (cq_s_polymw < 0.) {
2192 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2199 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2207 template<
typename TypeTag>
2208 std::optional<typename StandardWell<TypeTag>::Scalar>
2212 const SummaryState& summary_state,
2215 return computeBhpAtThpLimitProdWithAlq(simulator,
2217 this->getALQ(well_state),
2221 template<
typename TypeTag>
2222 std::optional<typename StandardWell<TypeTag>::Scalar>
2225 const SummaryState& summary_state,
2230 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2236 std::vector<Scalar> rates(3);
2237 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2238 this->adaptRatesForVFP(rates);
2242 Scalar max_pressure = 0.0;
2243 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2244 const int cell_idx = this->well_cells_[perf];
2245 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2246 const auto& fs = int_quants.fluidState();
2247 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2248 max_pressure = std::max(max_pressure, pressure_cell);
2253 this->connections_.rho(),
2255 this->getTHPConstraint(summary_state),
2259 auto v = frates(*bhpAtLimit);
2260 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2266 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2270 std::vector<Scalar> rates(3);
2271 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2272 this->adaptRatesForVFP(rates);
2279 this->connections_.rho(),
2281 this->getTHPConstraint(summary_state),
2287 auto v = frates(*bhpAtLimit);
2288 if (std::all_of(v.cbegin(), v.cend(), [](
Scalar i){ return i <= 0; }) ) {
2294 return std::nullopt;
2299 template<
typename TypeTag>
2300 std::optional<typename StandardWell<TypeTag>::Scalar>
2303 const SummaryState& summary_state,
2307 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2313 std::vector<Scalar> rates(3);
2314 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2320 this->connections_.rho(),
2331 template<
typename TypeTag>
2336 const Well::InjectionControls& inj_controls,
2337 const Well::ProductionControls& prod_controls,
2342 const int max_iter = this->param_.max_inner_iter_wells_;
2345 bool relax_convergence =
false;
2346 this->regularize_ =
false;
2348 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2350 if (it > this->param_.strict_inner_iter_wells_) {
2351 relax_convergence =
true;
2352 this->regularize_ =
true;
2355 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2357 converged = report.converged();
2363 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2370 initPrimaryVariablesEvaluation();
2371 }
while (it < max_iter);
2377 template<
typename TypeTag>
2382 const Well::InjectionControls& inj_controls,
2383 const Well::ProductionControls& prod_controls,
2387 const bool fixed_control ,
2388 const bool fixed_status )
2390 const int max_iter = this->param_.max_inner_iter_wells_;
2392 bool converged =
false;
2393 bool relax_convergence =
false;
2394 this->regularize_ =
false;
2395 const auto& summary_state = simulator.vanguard().summaryState();
2400 constexpr int min_its_after_switch = 4;
2401 int its_since_last_switch = min_its_after_switch;
2402 int switch_count= 0;
2404 const auto well_status_orig = this->wellStatus_;
2405 const auto operability_orig = this->operability_status_;
2406 auto well_status_cur = well_status_orig;
2407 int status_switch_count = 0;
2409 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2410 well_state.
well(this->index_of_well_).status == WellStatus::OPEN;
2412 const bool allow_switching =
2413 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2414 (!fixed_control || !fixed_status) && allow_open;
2416 bool changed =
false;
2417 bool final_check =
false;
2419 this->operability_status_.resetOperability();
2420 this->operability_status_.solvable =
true;
2422 its_since_last_switch++;
2423 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2424 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2425 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2426 inj_controls, prod_controls, wqTotal,
2427 deferred_logger, fixed_control, fixed_status);
2429 its_since_last_switch = 0;
2431 if (well_status_cur != this->wellStatus_) {
2432 well_status_cur = this->wellStatus_;
2433 status_switch_count++;
2436 if (!changed && final_check) {
2439 final_check =
false;
2443 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2445 if (it > this->param_.strict_inner_iter_wells_) {
2446 relax_convergence =
true;
2447 this->regularize_ =
true;
2450 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2452 converged = report.converged();
2456 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2458 its_since_last_switch = min_its_after_switch;
2465 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2466 initPrimaryVariablesEvaluation();
2468 }
while (it < max_iter);
2471 if (allow_switching){
2473 const bool is_stopped = this->wellIsStopped();
2474 if (this->wellHasTHPConstraints(summary_state)){
2475 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2476 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2478 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2482 this->wellStatus_ = well_status_orig;
2483 this->operability_status_ = operability_orig;
2484 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2485 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2486 deferred_logger.
debug(message);
2492 template<
typename TypeTag>
2493 std::vector<typename StandardWell<TypeTag>::Scalar>
2499 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2500 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2501 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2502 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2503 const int cell_idx = this->well_cells_[perf];
2504 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2505 std::vector<Scalar> mob(this->num_components_, 0.);
2506 getMobility(simulator, perf, mob, deferred_logger);
2507 std::vector<Scalar> cq_s(this->num_components_, 0.);
2508 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2509 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2510 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2512 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2513 cq_s, perf_rates, deferred_logger);
2514 for (
int comp = 0; comp < this->num_components_; ++comp) {
2515 well_q_s[comp] += cq_s[comp];
2518 const auto& comm = this->parallel_well_info_.communication();
2519 if (comm.size() > 1)
2521 comm.sum(well_q_s.data(), well_q_s.size());
2528 template <
typename TypeTag>
2529 std::vector<typename StandardWell<TypeTag>::Scalar>
2533 const int num_pri_vars = this->primary_variables_.numWellEq();
2534 std::vector<Scalar> retval(num_pri_vars);
2535 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2536 retval[ii] = this->primary_variables_.value(ii);
2545 template <
typename TypeTag>
2550 const int num_pri_vars = this->primary_variables_.numWellEq();
2551 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2552 this->primary_variables_.setValue(ii, it[ii]);
2554 return num_pri_vars;
2558 template <
typename TypeTag>
2562 const std::vector<EvalWell>& cq_s,
2563 const IntensiveQuantities& intQuants,
2566 auto fs = intQuants.fluidState();
2568 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2569 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2574 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2575 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2576 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2577 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2578 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2581 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2582 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2587 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2589 deferred_logger.
debug(
2590 fmt::format(
"Problematic d value {} obtained for well {}"
2591 " during calculateSinglePerf with rs {}"
2592 ", rv {}. Continue as if no dissolution (rs = 0) and"
2593 " vaporization (rv = 0) for this connection.",
2594 d, this->name(), fs.Rs(), fs.Rv()));
2595 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2597 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2598 cq_r_thermal = (cq_s[gasCompIdx] -
2599 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2600 (d * this->extendEval(fs.invB(phaseIdx)) );
2601 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2603 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2605 (d * this->extendEval(fs.invB(phaseIdx)) );
2611 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
2613 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2614 fs.setTemperature(this->well_ecl_.inj_temperature());
2615 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2616 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2617 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2618 paramCache.setRegionIndex(pvtRegionIdx);
2619 paramCache.setMaxOilSat(maxOilSaturation);
2620 paramCache.updatePhase(fs, phaseIdx);
2622 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2623 fs.setDensity(phaseIdx, rho);
2624 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2625 fs.setEnthalpy(phaseIdx, h);
2626 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2627 result += getValue(cq_r_thermal);
2630 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2631 result += Base::restrictEval(cq_r_thermal);
2635 return result * this->well_efficiency_factor_;
2639 template <
typename TypeTag>
2640 template<
class Value>
2642 StandardWell<TypeTag>::
2643 gasOilPerfRateInj(
const std::vector<Value>& cq_s,
2644 PerforationRates<Scalar>& perf_rates,
2647 const Value& pressure,
2649 DeferredLogger& deferred_logger)
const
2651 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2652 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2661 const Scalar d = 1.0 - getValue(rv) * getValue(rs);
2664 deferred_logger.debug(dValueError(d, this->name(),
2665 "gasOilPerfRateInj",
2670 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2673 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
2676 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2681 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2687 template <
typename TypeTag>
2688 template<
class Value>
2690 StandardWell<TypeTag>::
2691 gasOilPerfRateProd(std::vector<Value>& cq_s,
2692 PerforationRates<Scalar>& perf_rates,
2695 const Value& rvw)
const
2697 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2698 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2699 const Value cq_sOil = cq_s[oilCompIdx];
2700 const Value cq_sGas = cq_s[gasCompIdx];
2701 const Value dis_gas = rs * cq_sOil;
2702 const Value vap_oil = rv * cq_sGas;
2704 cq_s[gasCompIdx] += dis_gas;
2705 cq_s[oilCompIdx] += vap_oil;
2708 if (this->isProducer()) {
2709 perf_rates.dis_gas = getValue(dis_gas);
2710 perf_rates.vap_oil = getValue(vap_oil);
2713 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2714 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2715 const Value vap_wat = rvw * cq_sGas;
2716 cq_s[waterCompIdx] += vap_wat;
2717 if (this->isProducer())
2718 perf_rates.vap_wat = getValue(vap_wat);
2723 template <
typename TypeTag>
2724 template<
class Value>
2726 StandardWell<TypeTag>::
2727 gasWaterPerfRateProd(std::vector<Value>& cq_s,
2728 PerforationRates<Scalar>& perf_rates,
2730 const Value& rsw)
const
2732 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2733 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2734 const Value cq_sWat = cq_s[waterCompIdx];
2735 const Value cq_sGas = cq_s[gasCompIdx];
2736 const Value vap_wat = rvw * cq_sGas;
2737 const Value dis_gas_wat = rsw * cq_sWat;
2738 cq_s[waterCompIdx] += vap_wat;
2739 cq_s[gasCompIdx] += dis_gas_wat;
2740 if (this->isProducer()) {
2741 perf_rates.vap_wat = getValue(vap_wat);
2742 perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
2747 template <
typename TypeTag>
2748 template<
class Value>
2750 StandardWell<TypeTag>::
2751 gasWaterPerfRateInj(
const std::vector<Value>& cq_s,
2752 PerforationRates<Scalar>& perf_rates,
2755 const Value& pressure,
2756 DeferredLogger& deferred_logger)
const
2759 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2760 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2762 const Scalar dw = 1.0 - getValue(rvw) * getValue(rsw);
2765 deferred_logger.debug(dValueError(dw, this->name(),
2766 "gasWaterPerfRateInj",
2767 rsw, rvw, pressure));
2771 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
2774 perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
2779 template <
typename TypeTag>
2780 template<
class Value>
2782 StandardWell<TypeTag>::
2783 disOilVapWatVolumeRatio(Value& volumeRatio,
2786 const Value& pressure,
2787 const std::vector<Value>& cmix_s,
2788 const std::vector<Value>& b_perfcells_dense,
2789 DeferredLogger& deferred_logger)
const
2791 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2792 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2794 const Value d = 1.0 - rvw * rsw;
2797 deferred_logger.debug(dValueError(d, this->name(),
2798 "disOilVapWatVolumeRatio",
2799 rsw, rvw, pressure));
2801 const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
2802 : cmix_s[waterCompIdx];
2803 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
2805 const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
2806 : cmix_s[gasCompIdx];
2807 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
2811 template <
typename TypeTag>
2812 template<
class Value>
2814 StandardWell<TypeTag>::
2815 gasOilVolumeRatio(Value& volumeRatio,
2818 const Value& pressure,
2819 const std::vector<Value>& cmix_s,
2820 const std::vector<Value>& b_perfcells_dense,
2821 DeferredLogger& deferred_logger)
const
2823 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2824 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2826 const Value d = 1.0 - rv * rs;
2829 deferred_logger.debug(dValueError(d, this->name(),
2830 "gasOilVolumeRatio",
2833 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
2834 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
2836 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
2837 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:38
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:186
Class handling assemble of the equation system for StandardWell.
Definition: StandardWellAssemble.hpp:43
Scalar pressure_diff(const unsigned perf) const
Returns pressure drop for a given perforation.
Definition: StandardWellConnections.hpp:105
StdWellConnections connections_
Connection level values.
Definition: StandardWellEval.hpp:107
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:1365
EvalWell wpolymermw(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2017
void computeWellPotentials(const Simulator &simulator, const WellState< Scalar > &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) override
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1742
typename StdWellEval::EvalWell EvalWell
Definition: StandardWell.hpp:119
std::optional< Scalar > computeBhpAtThpLimitProdWithAlq(const Simulator &ebos_simulator, const SummaryState &summary_state, const Scalar alq_value, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2224
typename StdWellEval::BVectorWell BVectorWell
Definition: StandardWell.hpp:120
void checkOperabilityUnderBHPLimit(const WellState< Scalar > &well_state, const Simulator &simulator, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:937
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:2334
void computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1725
void addWellContributions(SparseMatrixAdapter &mat) const override
Definition: StandardWell_impl.hpp:1912
std::vector< Scalar > getPrimaryVars() const override
Definition: StandardWell_impl.hpp:2531
void initPrimaryVariablesEvaluation() override
Definition: StandardWell_impl.hpp:113
void updateIPRImplicit(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:861
void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator< Scalar > &wellPICalc, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1239
void updateWaterMobilityWithPolymer(const Simulator &simulator, const int perf, std::vector< EvalWell > &mob_water, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1848
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:2380
StandardWell(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters ¶m, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Definition: StandardWell_impl.hpp:72
void updateWaterThroughput(const double dt, WellState< Scalar > &well_state) const override
Definition: StandardWell_impl.hpp:2046
std::vector< Scalar > computeCurrentWellRates(const Simulator &ebosSimulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2495
virtual void init(const PhaseUsage *phase_usage_arg, const std::vector< Scalar > &depth_arg, const Scalar 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 updatePrimaryVariablesNewton(const BVectorWell &dwells, const bool stop_or_zero_rate_target, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:725
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:79
void calculateExplicitQuantities(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1405
void solveEqAndUpdateWellState(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1383
bool computeWellPotentialsImplicit(const Simulator &ebos_simulator, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1618
void recoverWellSolutionAndUpdateWellState(const Simulator &simulator, const BVector &x, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1452
typename StdWellEval::StdWellConnections::Properties WellConnectionProps
Definition: StandardWell.hpp:270
void computeWellRatesWithBhp(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1472
WellConnectionProps computePropertiesForWellConnectionPressures(const Simulator &simulator, const WellState< Scalar > &well_state) const
Definition: StandardWell_impl.hpp:1144
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:647
void computeWellRatesWithBhpIterations(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1518
void updateIPR(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:768
std::vector< Scalar > computeWellPotentialWithTHP(const Simulator &ebosSimulator, DeferredLogger &deferred_logger, const WellState< Scalar > &well_state) const
Definition: StandardWell_impl.hpp:1583
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:126
void updateWellState(const Simulator &simulator, const BVectorWell &dwells, WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:703
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:344
void computeWellConnectionDensitesPressures(const Simulator &simulator, const WellState< Scalar > &well_state, const WellConnectionProps &props, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1314
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1420
std::optional< Scalar > computeBhpAtThpLimitProd(const WellState< Scalar > &well_state, const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2210
void checkConvergenceExtraEqs(const std::vector< Scalar > &res, ConvergenceReport &report) const
Definition: StandardWell_impl.hpp:2148
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:1696
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:2167
bool openCrossFlowAvoidSingularity(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1133
void updatePrimaryVariables(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1814
virtual ConvergenceReport getWellConvergence(const Simulator &simulator, const WellState< Scalar > &well_state, const std::vector< Scalar > &B_avg, DeferredLogger &deferred_logger, const bool relax_tolerance) const override
check whether the well equations get converged for this well
Definition: StandardWell_impl.hpp:1195
std::optional< Scalar > computeBhpAtThpLimitInj(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2302
bool allDrawDownWrongDirection(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1055
EvalWell pskin(const Scalar throughput, const EvalWell &water_velocity, const EvalWell &poly_inj_conc, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1972
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:748
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:2100
int setPrimaryVars(typename std::vector< Scalar >::const_iterator it) override
Definition: StandardWell_impl.hpp:2548
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:370
void checkOperabilityUnderTHPLimit(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1006
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:1920
EvalWell pskinwater(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1940
bool canProduceInjectWithCurrentBhp(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1097
void handleInjectivityRate(const Simulator &simulator, const int perf, std::vector< EvalWell > &cq_s) const
Definition: StandardWell_impl.hpp:2077
Scalar getRefDensity() const override
Definition: StandardWell_impl.hpp:1837
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:497
Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const override
Definition: StandardWell_impl.hpp:1799
EvalWell getQs(const int compIdx) const
Returns scaled rate for a component.
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
std::optional< Scalar > computeBhpAtThpLimitProd(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar maxPerfPress, const Scalar rho, const Scalar alq_value, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for a producer.
std::optional< Scalar > computeBhpAtThpLimitInj(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar rho, const Scalar flo_rel_tol, const int max_iteration, const bool throwOnError, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for an injector.
Scalar calculateThpFromBhp(const std::vector< Scalar > &rates, const Scalar bhp, const Scalar rho, const std::optional< Scalar > &alq, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Calculates THP from BHP.
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
Definition: WellConvergence.hpp:38
void prepareForPotentialCalculations(const SummaryState &summary_state, WellState< Scalar > &well_state, Well::InjectionControls &inj_controls, Well::ProductionControls &prod_controls) const
std::pair< bool, bool > computeWellPotentials(std::vector< Scalar > &well_potentials, const WellState< Scalar > &well_state)
const int num_components_
Definition: WellInterfaceGeneric.hpp:282
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:273
Definition: WellInterfaceIndices.hpp:34
Definition: WellInterface.hpp:73
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:78
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1773
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:100
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, 1, 1 > > PressureMatrix
Definition: WellInterface.hpp:97
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:82
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:79
Dune::BlockVector< VectorBlockType > BVector
Definition: WellInterface.hpp:96
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:80
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:81
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:514
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: WellInterface.hpp:84
Definition: WellProdIndexCalculator.hpp:37
Scalar connectionProdIndStandard(const std::size_t connIdx, const Scalar connMobility) const
Definition: WellState.hpp:62
const SingleWellState< Scalar > & well(std::size_t well_index) const
Definition: WellState.hpp:307
std::vector< Scalar > & wellRates(std::size_t well_index)
One rate per well and phase.
Definition: WellState.hpp:272
@ NONE
Definition: DeferredLogger.hpp:46
Definition: blackoilboundaryratevector.hh:37
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:147
Static data associated with a well perforation.
Definition: PerforationData.hpp:30
Definition: PerforationData.hpp:41
Scalar dis_gas
Definition: PerforationData.hpp:42
Scalar vap_wat
Definition: PerforationData.hpp:45
Scalar vap_oil
Definition: PerforationData.hpp:44
Scalar dis_gas_in_water
Definition: PerforationData.hpp:43
Definition: BlackoilPhases.hpp:46