22#include <opm/common/Exceptions.hpp>
24#include <opm/input/eclipse/Units/Units.hpp>
26#include <opm/material/densead/EvaluationFormat.hpp>
28#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
29#include <opm/simulators/wells/StandardWellAssemble.hpp>
30#include <opm/simulators/wells/VFPHelpers.hpp>
31#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
32#include <opm/simulators/wells/WellConvergence.hpp>
34#include <fmt/format.h>
43template<
class dValue,
class Value>
44auto dValueError(
const dValue& d,
45 const std::string& name,
46 const std::string& methodName,
49 const Value& pressure)
51 return fmt::format(
"Problematic d value {} obtained for well {}"
52 " during {} calculations with rs {}"
53 ", rv {} and pressure {}."
54 " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
55 " for this connection.", d, name, methodName, Rs, Rv, pressure);
63 template<
typename TypeTag>
64 StandardWell<TypeTag>::
65 StandardWell(
const Well& well,
66 const ParallelWellInfo& pw_info,
68 const ModelParameters& param,
69 const RateConverterType& rate_converter,
70 const int pvtRegionIdx,
71 const int num_components,
73 const int index_of_well,
74 const std::vector<PerforationData>& perf_data)
75 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
76 , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
79 assert(this->num_components_ == numWellConservationEq);
86 template<
typename TypeTag>
88 StandardWell<TypeTag>::
89 init(
const PhaseUsage* phase_usage_arg,
90 const std::vector<double>& depth_arg,
91 const double gravity_arg,
93 const std::vector< Scalar >& B_avg,
94 const bool changed_to_open_this_step)
96 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
97 this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
104 template<
typename TypeTag>
105 void StandardWell<TypeTag>::
106 initPrimaryVariablesEvaluation()
108 this->primary_variables_.init();
115 template<
typename TypeTag>
116 template<
class Value>
118 StandardWell<TypeTag>::
119 computePerfRate(
const IntensiveQuantities& intQuants,
120 const std::vector<Value>& mob,
125 std::vector<Value>& cq_s,
126 PerforationRates& perf_rates,
127 DeferredLogger& deferred_logger)
const
129 auto obtain = [
this](
const Eval& value)
131 if constexpr (std::is_same_v<Value, Scalar>) {
132 static_cast<void>(
this);
133 return getValue(value);
135 return this->extendEval(value);
138 auto obtainN = [](
const auto& value)
140 if constexpr (std::is_same_v<Value, Scalar>) {
141 return getValue(value);
146 auto zeroElem = [
this]()
148 if constexpr (std::is_same_v<Value, Scalar>) {
149 static_cast<void>(
this);
152 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
156 const auto& fs = intQuants.fluidState();
157 const Value pressure = obtain(this->getPerfCellPressure(fs));
158 const Value rs = obtain(fs.Rs());
159 const Value rv = obtain(fs.Rv());
160 const Value rvw = obtain(fs.Rvw());
161 const Value rsw = obtain(fs.Rsw());
163 std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
164 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
165 if (!FluidSystem::phaseIsActive(phaseIdx)) {
168 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
169 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
171 if constexpr (has_solvent) {
172 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
175 if constexpr (has_zFraction) {
176 if (this->isInjector()) {
177 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
178 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
179 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
183 Value skin_pressure = zeroElem();
185 if (this->isInjector()) {
186 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
187 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
192 std::vector<Value> cmix_s(this->numComponents(), zeroElem());
193 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
194 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
215 template<
typename TypeTag>
216 template<
class Value>
218 StandardWell<TypeTag>::
219 computePerfRate(
const std::vector<Value>& mob,
220 const Value& pressure,
226 std::vector<Value>& b_perfcells_dense,
230 const Value& skin_pressure,
231 const std::vector<Value>& cmix_s,
232 std::vector<Value>& cq_s,
233 PerforationRates& perf_rates,
234 DeferredLogger& deferred_logger)
const
237 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
238 Value drawdown = pressure - well_pressure;
239 if (this->isInjector()) {
240 drawdown += skin_pressure;
246 if (!allow_cf && this->isInjector()) {
251 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
252 const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
253 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
256 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
257 gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw);
258 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
259 gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw);
263 if (!allow_cf && this->isProducer()) {
268 Value total_mob_dense = mob[0];
269 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
270 total_mob_dense += mob[componentIdx];
274 const Value cqt_i = - Tw * (total_mob_dense * drawdown);
277 Value volumeRatio = bhp * 0.0;
279 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
280 disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
281 cmix_s, b_perfcells_dense, deferred_logger);
283 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
284 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
285 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
289 if constexpr (Indices::enableSolvent) {
290 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
293 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
294 gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
295 cmix_s, b_perfcells_dense, deferred_logger);
298 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
299 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
300 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
302 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
303 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
304 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
309 Value cqt_is = cqt_i / volumeRatio;
310 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
311 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
315 if (this->isProducer()) {
316 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
317 gasOilPerfRateInj(cq_s, perf_rates,
318 rv, rs, pressure, rvw, deferred_logger);
320 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
322 gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
323 pressure, deferred_logger);
330 template<
typename TypeTag>
332 StandardWell<TypeTag>::
333 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
335 const Well::InjectionControls& inj_controls,
336 const Well::ProductionControls& prod_controls,
337 WellState& well_state,
338 const GroupState& group_state,
339 DeferredLogger& deferred_logger)
343 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
346 this->linSys_.clear();
348 assembleWellEqWithoutIterationImpl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
354 template<
typename TypeTag>
356 StandardWell<TypeTag>::
357 assembleWellEqWithoutIterationImpl(
const Simulator& ebosSimulator,
359 const Well::InjectionControls& inj_controls,
360 const Well::ProductionControls& prod_controls,
361 WellState& well_state,
362 const GroupState& group_state,
363 DeferredLogger& deferred_logger)
366 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
367 const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
369 auto& ws = well_state.well(this->index_of_well_);
370 ws.phase_mixing_rates.fill(0.0);
372 const int np = this->number_of_phases_;
374 std::vector<RateVector> connectionRates = this->connectionRates_;
375 auto& perf_data = ws.perf_data;
376 auto& perf_rates = perf_data.phase_rates;
377 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
379 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
380 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
381 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
382 calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
385 if constexpr (has_polymer && Base::has_polymermw) {
386 if (this->isInjector()) {
387 handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
390 const int cell_idx = this->well_cells_[perf];
391 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
393 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
395 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
397 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
398 assemblePerforationEq(cq_s_effective,
401 this->primary_variables_.numWellEq(),
405 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
406 auto& perf_rate_solvent = perf_data.solvent_rates;
407 perf_rate_solvent[perf] = cq_s[componentIdx].value();
409 perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
413 if constexpr (has_zFraction) {
414 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
415 assembleZFracEq(cq_s_zfrac_effective,
417 this->primary_variables_.numWellEq(),
422 this->connectionRates_ = connectionRates;
427 const auto& comm = this->parallel_well_info_.communication();
428 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
432 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
435 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
438 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
439 if (FluidSystem::numActivePhases() > 1) {
441 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
442 this->F0_[componentIdx]) * volume / dt;
444 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
445 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
446 assembleSourceEq(resWell_loc,
448 this->primary_variables_.numWellEq(),
452 const auto& summaryState = ebosSimulator.vanguard().summaryState();
453 const Schedule& schedule = ebosSimulator.vanguard().schedule();
454 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
455 assembleControlEq(well_state, group_state,
456 schedule, summaryState,
457 inj_controls, prod_controls,
458 this->primary_variables_,
459 this->connections_.rho(),
466 this->linSys_.invert();
468 OPM_DEFLOG_THROW(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
475 template<
typename TypeTag>
477 StandardWell<TypeTag>::
478 calculateSinglePerf(
const Simulator& ebosSimulator,
480 WellState& well_state,
481 std::vector<RateVector>& connectionRates,
482 std::vector<EvalWell>& cq_s,
483 EvalWell& water_flux_s,
484 EvalWell& cq_s_zfrac_effective,
485 DeferredLogger& deferred_logger)
const
487 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
488 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
489 const int cell_idx = this->well_cells_[perf];
490 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
491 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
492 getMobility(ebosSimulator, perf, mob, deferred_logger);
494 PerforationRates perf_rates;
495 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
496 const double Tw = this->well_index_[perf] * trans_mult;
497 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
498 cq_s, perf_rates, deferred_logger);
500 auto& ws = well_state.well(this->index_of_well_);
501 auto& perf_data = ws.perf_data;
502 if constexpr (has_polymer && Base::has_polymermw) {
503 if (this->isInjector()) {
506 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
507 water_flux_s = cq_s[water_comp_idx];
510 handleInjectivityRate(ebosSimulator, perf, cq_s);
515 if (this->isProducer()) {
516 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
517 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
518 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
519 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
522 if constexpr (has_energy) {
523 connectionRates[perf][Indices::contiEnergyEqIdx] =
524 connectionRateEnergy(ebosSimulator.problem().maxOilSaturation(cell_idx),
525 cq_s, intQuants, deferred_logger);
528 if constexpr (has_polymer) {
529 std::variant<Scalar,EvalWell> polymerConcentration;
530 if (this->isInjector()) {
531 polymerConcentration = this->wpolymer();
533 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
534 intQuants.polymerViscosityCorrection());
537 [[maybe_unused]] EvalWell cq_s_poly;
538 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
540 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
541 cq_s, polymerConcentration);
543 if constexpr (Base::has_polymermw) {
544 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
545 perf, connectionRates, deferred_logger);
549 if constexpr (has_foam) {
550 std::variant<Scalar,EvalWell> foamConcentration;
551 if (this->isInjector()) {
552 foamConcentration = this->wfoam();
554 foamConcentration = this->extendEval(intQuants.foamConcentration());
556 connectionRates[perf][Indices::contiFoamEqIdx] =
557 this->connections_.connectionRateFoam(cq_s, foamConcentration,
558 FoamModule::transportPhase(),
562 if constexpr (has_zFraction) {
563 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
564 if (this->isInjector()) {
565 solventConcentration = this->wsolvent();
567 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
568 this->extendEval(intQuants.yVolume())};
570 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
571 cq_s_zfrac_effective) =
572 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
573 perf_rates.dis_gas, cq_s,
574 solventConcentration);
577 if constexpr (has_brine) {
578 std::variant<Scalar,EvalWell> saltConcentration;
579 if (this->isInjector()) {
580 saltConcentration = this->wsalt();
582 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
585 connectionRates[perf][Indices::contiBrineEqIdx] =
586 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
587 perf_rates.vap_wat, cq_s,
591 if constexpr (has_micp) {
592 std::variant<Scalar,EvalWell> microbialConcentration;
593 std::variant<Scalar,EvalWell> oxygenConcentration;
594 std::variant<Scalar,EvalWell> ureaConcentration;
595 if (this->isInjector()) {
596 microbialConcentration = this->wmicrobes();
597 oxygenConcentration = this->woxygen();
598 ureaConcentration = this->wurea();
600 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
601 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
602 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
604 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
605 connectionRates[perf][Indices::contiOxygenEqIdx],
606 connectionRates[perf][Indices::contiUreaEqIdx]) =
607 this->connections_.connectionRatesMICP(cq_s,
608 microbialConcentration,
614 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
619 template<
typename TypeTag>
620 template<
class Value>
622 StandardWell<TypeTag>::
623 getMobility(
const Simulator& ebosSimulator,
625 std::vector<Value>& mob,
626 DeferredLogger& deferred_logger)
const
628 auto obtain = [
this](
const Eval& value)
630 if constexpr (std::is_same_v<Value, Scalar>) {
631 static_cast<void>(
this);
632 return getValue(value);
634 return this->extendEval(value);
637 WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob,
638 obtain, deferred_logger);
641 if constexpr (has_polymer) {
642 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
643 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
648 if constexpr (!Base::has_polymermw) {
649 if constexpr (std::is_same_v<Value, Scalar>) {
650 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
651 for (std::size_t i = 0; i < mob.size(); ++i) {
652 mob_eval[i].setValue(mob[i]);
654 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
655 for (std::size_t i = 0; i < mob.size(); ++i) {
656 mob[i] = getValue(mob_eval[i]);
659 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
665 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
666 const double bhp = this->primary_variables_.value(Bhp);
667 const double perf_press = bhp + this->connections_.pressure_diff(perf);
668 const double multiplier = this->getInjMult(perf, bhp, perf_press);
669 for (std::size_t i = 0; i < mob.size(); ++i) {
670 mob[i] *= multiplier;
676 template<
typename TypeTag>
678 StandardWell<TypeTag>::
679 updateWellState(
const SummaryState& summary_state,
680 const BVectorWell& dwells,
681 WellState& well_state,
682 DeferredLogger& deferred_logger)
684 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
686 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
687 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
689 updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
690 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
697 template<
typename TypeTag>
699 StandardWell<TypeTag>::
700 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
701 const bool stop_or_zero_rate_target,
702 DeferredLogger& deferred_logger)
704 const double dFLimit = this->param_.dwell_fraction_max_;
705 const double dBHPLimit = this->param_.dbhp_max_rel_;
706 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit);
709 if constexpr (Base::has_polymermw) {
710 this->primary_variables_.updateNewtonPolyMW(dwells);
713 this->primary_variables_.checkFinite(deferred_logger);
720 template<
typename TypeTag>
722 StandardWell<TypeTag>::
723 updateWellStateFromPrimaryVariables(
const bool stop_or_zero_rate_target,
724 WellState& well_state,
725 const SummaryState& summary_state,
726 DeferredLogger& deferred_logger)
const
728 this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
731 if constexpr (Base::has_polymermw) {
732 this->primary_variables_.copyToWellStatePolyMW(well_state);
740 template<
typename TypeTag>
742 StandardWell<TypeTag>::
743 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
748 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
749 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
751 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
752 std::vector<Scalar> mob(this->num_components_, 0.0);
753 getMobility(ebos_simulator, perf, mob, deferred_logger);
755 const int cell_idx = this->well_cells_[perf];
756 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
757 const auto& fs = int_quantities.fluidState();
759 double p_r = this->getPerfCellPressure(fs).value();
762 std::vector<double> b_perf(this->num_components_);
763 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
764 if (!FluidSystem::phaseIsActive(phase)) {
767 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
768 b_perf[comp_idx] = fs.invB(phase).value();
770 if constexpr (has_solvent) {
771 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
775 const double h_perf = this->connections_.pressure_diff(perf);
776 const double pressure_diff = p_r - h_perf;
781 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
782 deferred_logger.debug(
"CROSSFLOW_IPR",
783 "cross flow found when updateIPR for well " + name()
784 +
" . The connection is ignored in IPR calculations");
790 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
792 std::vector<double> ipr_a_perf(this->ipr_a_.size());
793 std::vector<double> ipr_b_perf(this->ipr_b_.size());
794 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
795 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
796 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
797 ipr_b_perf[comp_idx] += tw_mob;
801 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
802 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
803 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
804 const double rs = (fs.Rs()).value();
805 const double rv = (fs.Rv()).value();
807 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
808 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
810 ipr_a_perf[gas_comp_idx] += dis_gas_a;
811 ipr_a_perf[oil_comp_idx] += vap_oil_a;
813 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
814 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
816 ipr_b_perf[gas_comp_idx] += dis_gas_b;
817 ipr_b_perf[oil_comp_idx] += vap_oil_b;
820 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
821 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
822 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
825 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
826 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
830 template<
typename TypeTag>
832 StandardWell<TypeTag>::
833 checkOperabilityUnderBHPLimit(
const WellState& well_state,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
835 const auto& summaryState = ebos_simulator.vanguard().summaryState();
836 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
839 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
840 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
843 double total_ipr_mass_rate = 0.0;
844 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
846 if (!FluidSystem::phaseIsActive(phaseIdx)) {
850 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
851 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
853 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
854 total_ipr_mass_rate += ipr_rate * rho;
856 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
857 this->operability_status_.operable_under_only_bhp_limit =
false;
861 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
865 std::vector<double> well_rates_bhp_limit;
866 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
868 this->adaptRatesForVFP(well_rates_bhp_limit);
869 const double thp_limit = this->getTHPConstraint(summaryState);
870 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
872 this->connections_.rho(),
873 this->getALQ(well_state),
876 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
877 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
888 this->operability_status_.operable_under_only_bhp_limit =
true;
889 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
897 template<
typename TypeTag>
899 StandardWell<TypeTag>::
900 checkOperabilityUnderTHPLimit(
const Simulator& ebos_simulator,
const WellState& well_state, DeferredLogger& deferred_logger)
902 const auto& summaryState = ebos_simulator.vanguard().summaryState();
903 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger)
904 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
907 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
909 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
910 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
911 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
913 const double thp_limit = this->getTHPConstraint(summaryState);
914 if (this->isProducer() && *obtain_bhp < thp_limit) {
915 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
916 +
" bars is SMALLER than thp limit "
917 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
918 +
" bars as a producer for well " + name();
919 deferred_logger.debug(msg);
921 else if (this->isInjector() && *obtain_bhp > thp_limit) {
922 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
923 +
" bars is LARGER than thp limit "
924 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
925 +
" bars as a injector for well " + name();
926 deferred_logger.debug(msg);
929 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
930 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
931 if (!this->wellIsStopped()) {
932 const double thp_limit = this->getTHPConstraint(summaryState);
933 deferred_logger.debug(
" could not find bhp value at thp limit "
934 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
935 +
" bar for well " + name() +
", the well might need to be closed ");
944 template<
typename TypeTag>
946 StandardWell<TypeTag>::
947 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
949 bool all_drawdown_wrong_direction =
true;
951 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
952 const int cell_idx = this->well_cells_[perf];
953 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
954 const auto& fs = intQuants.fluidState();
956 const double pressure = this->getPerfCellPressure(fs).value();
957 const double bhp = this->primary_variables_.eval(Bhp).value();
960 const double well_pressure = bhp + this->connections_.pressure_diff(perf);
961 const double drawdown = pressure - well_pressure;
966 if ( (drawdown < 0. && this->isInjector()) ||
967 (drawdown > 0. && this->isProducer()) ) {
968 all_drawdown_wrong_direction =
false;
973 const auto& comm = this->parallel_well_info_.communication();
976 all_drawdown_wrong_direction =
977 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
980 return all_drawdown_wrong_direction;
986 template<
typename TypeTag>
988 StandardWell<TypeTag>::
989 canProduceInjectWithCurrentBhp(
const Simulator& ebos_simulator,
990 const WellState& well_state,
991 DeferredLogger& deferred_logger)
993 const double bhp = well_state.well(this->index_of_well_).bhp;
994 std::vector<double> well_rates;
995 computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
997 const double sign = (this->isProducer()) ? -1. : 1.;
998 const double threshold = sign * std::numeric_limits<double>::min();
1000 bool can_produce_inject =
false;
1001 for (
const auto value : well_rates) {
1002 if (this->isProducer() && value < threshold) {
1003 can_produce_inject =
true;
1005 }
else if (this->isInjector() && value > threshold) {
1006 can_produce_inject =
true;
1011 if (!can_produce_inject) {
1012 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1015 return can_produce_inject;
1022 template<
typename TypeTag>
1024 StandardWell<TypeTag>::
1025 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1027 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1033 template<
typename TypeTag>
1035 StandardWell<TypeTag>::
1036 computePropertiesForWellConnectionPressures(
const Simulator& ebosSimulator,
1037 const WellState& well_state,
1038 WellConnectionProps& props)
const
1040 std::function<Scalar(
int,
int)> getTemperature =
1041 [&ebosSimulator](
int cell_idx,
int phase_idx)
1043 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().temperature(phase_idx).value();
1045 std::function<Scalar(
int)> getSaltConcentration =
1046 [&ebosSimulator](
int cell_idx)
1048 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().saltConcentration().value();
1050 std::function<int(
int)> getPvtRegionIdx =
1051 [&ebosSimulator](
int cell_idx)
1053 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().pvtRegionIndex();
1055 std::function<Scalar(
int)> getInvFac =
1056 [&ebosSimulator](
int cell_idx)
1058 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1060 std::function<Scalar(
int)> getSolventDensity =
1061 [&ebosSimulator](
int cell_idx)
1063 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventRefDensity();
1066 this->connections_.computePropertiesForPressures(well_state,
1068 getSaltConcentration,
1079 template<
typename TypeTag>
1084 const std::vector<double>&
B_avg,
1090 assert((
int(
B_avg.size()) ==
this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1095 this->param_.tolerance_wells_ *
stricter_factor : this->param_.tolerance_wells_;
1097 std::vector<double>
res;
1100 this->param_.max_residual_allowed_,
1102 this->param_.relaxed_tolerance_flow_well_,
1107 checkConvergenceExtraEqs(
res, report);
1116 template<
typename TypeTag>
1127 return ebosSimulator.model()
1128 .intensiveQuantities(
cell_idx, 0).fluidState();
1131 const int np = this->number_of_phases_;
1132 auto setToZero = [np](
double* x) ->
void
1134 std::fill_n(x, np, 0.0);
1137 auto addVector = [np](
const double* src,
double* dest) ->
void
1139 std::transform(src, src + np, dest, dest, std::plus<>{});
1142 auto& ws = well_state.well(this->index_of_well_);
1143 auto& perf_data = ws.perf_data;
1144 auto* wellPI = ws.productivity_index.data();
1145 auto* connPI = perf_data.prod_index.data();
1149 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1150 auto subsetPerfID = 0;
1152 for (
const auto& perf : *this->perf_data_) {
1153 auto allPerfID = perf.ecl_index;
1155 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1160 std::vector<Scalar> mob(this->num_components_, 0.0);
1161 getMobility(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1163 const auto& fs = fluidState(subsetPerfID);
1166 if (this->isInjector()) {
1167 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1168 mob, connPI, deferred_logger);
1171 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1174 addVector(connPI, wellPI);
1181 const auto& comm = this->parallel_well_info_.communication();
1182 if (comm.size() > 1) {
1183 comm.sum(wellPI, np);
1186 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1187 "Internal logic error in processing connections for PI/II");
1192 template<
typename TypeTag>
1194 StandardWell<TypeTag>::
1195 computeWellConnectionDensitesPressures(
const Simulator& ebosSimulator,
1196 const WellState& well_state,
1197 const WellConnectionProps& props,
1198 DeferredLogger& deferred_logger)
1200 std::function<Scalar(
int,
int)> invB =
1201 [&ebosSimulator](
int cell_idx,
int phase_idx)
1203 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().invB(phase_idx).value();
1205 std::function<Scalar(
int,
int)> mobility =
1206 [&ebosSimulator](
int cell_idx,
int phase_idx)
1208 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).mobility(phase_idx).value();
1210 std::function<Scalar(
int)> invFac =
1211 [&ebosSimulator](
int cell_idx)
1213 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1215 std::function<Scalar(
int)> solventMobility =
1216 [&ebosSimulator](
int cell_idx)
1218 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventMobility().value();
1221 this->connections_.computeProperties(well_state,
1234 template<
typename TypeTag>
1236 StandardWell<TypeTag>::
1237 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1238 const WellState& well_state,
1239 DeferredLogger& deferred_logger)
1244 WellConnectionProps props;
1245 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, props);
1246 computeWellConnectionDensitesPressures(ebosSimulator, well_state,
1247 props, deferred_logger);
1254 template<
typename TypeTag>
1256 StandardWell<TypeTag>::
1257 solveEqAndUpdateWellState(
const SummaryState& summary_state,
1258 WellState& well_state,
1259 DeferredLogger& deferred_logger)
1261 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1265 BVectorWell dx_well(1);
1266 dx_well[0].resize(this->primary_variables_.numWellEq());
1267 this->linSys_.solve( dx_well);
1269 updateWellState(summary_state, dx_well, well_state, deferred_logger);
1276 template<
typename TypeTag>
1278 StandardWell<TypeTag>::
1279 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1280 const WellState& well_state,
1281 DeferredLogger& deferred_logger)
1283 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1284 updatePrimaryVariables(summary_state, well_state, deferred_logger);
1285 initPrimaryVariablesEvaluation();
1286 computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
1287 this->computeAccumWell();
1292 template<
typename TypeTag>
1295 apply(
const BVector& x, BVector&
Ax)
const
1297 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1299 if (this->param_.matrix_add_well_contributions_)
1311 template<
typename TypeTag>
1316 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1324 template<
typename TypeTag>
1332 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1335 xw[0].resize(this->primary_variables_.numWellEq());
1337 this->linSys_.recoverSolutionWell(x,
xw);
1344 template<
typename TypeTag>
1353 const int np = this->number_of_phases_;
1356 const bool allow_cf = this->getAllowCrossFlow();
1358 for (
int perf = 0;
perf < this->number_of_perforations_; ++
perf) {
1362 std::vector<Scalar>
mob(this->num_components_, 0.);
1367 std::vector<Scalar>
cq_s(this->num_components_, 0.);
1372 for(
int p = 0;
p <
np; ++
p) {
1377 if constexpr (has_solvent) {
1378 const auto& pu = this->phaseUsage();
1379 assert(pu.phase_used[Gas]);
1380 const int gas_pos = pu.phase_pos[Gas];
1381 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1384 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1389 template<
typename TypeTag>
1391 StandardWell<TypeTag>::
1392 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
1394 std::vector<double>& well_flux,
1395 DeferredLogger& deferred_logger)
const
1399 StandardWell<TypeTag> well_copy(*
this);
1404 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1405 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1408 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1409 auto inj_controls = well_copy.well_ecl_.isInjector()
1410 ? well_copy.well_ecl_.injectionControls(summary_state)
1411 : Well::InjectionControls(0);
1412 auto prod_controls = well_copy.well_ecl_.isProducer()
1413 ? well_copy.well_ecl_.productionControls(summary_state) :
1414 Well::ProductionControls(0);
1417 auto& ws = well_state_copy.well(this->index_of_well_);
1418 if (well_copy.well_ecl_.isInjector()) {
1419 inj_controls.bhp_limit = bhp;
1420 ws.injection_cmode = Well::InjectorCMode::BHP;
1422 prod_controls.bhp_limit = bhp;
1423 ws.production_cmode = Well::ProducerCMode::BHP;
1428 const int np = this->number_of_phases_;
1429 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1430 for (
int phase = 0; phase < np; ++phase){
1431 well_state_copy.wellRates(this->index_of_well_)[phase]
1432 = sign * ws.well_potentials[phase];
1434 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1436 const double dt = ebosSimulator.timeStepSize();
1437 const bool converged = well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1439 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1440 " potentials are computed based on unconverged solution";
1441 deferred_logger.debug(msg);
1443 well_copy.updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
1444 well_copy.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
1445 well_copy.initPrimaryVariablesEvaluation();
1446 well_copy.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1452 template<
typename TypeTag>
1454 StandardWell<TypeTag>::
1455 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
1456 DeferredLogger& deferred_logger,
1457 const WellState &well_state)
const
1459 std::vector<double> potentials(this->number_of_phases_, 0.0);
1460 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1462 const auto& well = this->well_ecl_;
1463 if (well.isInjector()){
1464 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1465 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1466 if (bhp_at_thp_limit) {
1467 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1468 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1470 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1471 "Failed in getting converged thp based potential calculation for well "
1472 + name() +
". Instead the bhp based value is used");
1473 const double bhp = controls.bhp_limit;
1474 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1477 computeWellRatesWithThpAlqProd(
1478 ebos_simulator, summary_state,
1479 deferred_logger, potentials, this->getALQ(well_state)
1488 template<
typename TypeTag>
1490 StandardWell<TypeTag>::
1491 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &ebos_simulator,
1492 const SummaryState &summary_state,
1493 DeferredLogger &deferred_logger,
1494 std::vector<double> &potentials,
1498 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1499 ebos_simulator, summary_state, alq, deferred_logger);
1500 if (bhp_at_thp_limit) {
1501 const auto& controls = this->well_ecl_.productionControls(summary_state);
1502 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1503 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1506 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1507 "Failed in getting converged thp based potential calculation for well "
1508 + name() +
". Instead the bhp based value is used");
1509 const auto& controls = this->well_ecl_.productionControls(summary_state);
1510 bhp = controls.bhp_limit;
1511 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1516 template<
typename TypeTag>
1518 StandardWell<TypeTag>::
1519 computeWellRatesWithThpAlqProd(
const Simulator &ebos_simulator,
1520 const SummaryState &summary_state,
1521 DeferredLogger &deferred_logger,
1522 std::vector<double> &potentials,
1526 computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1533 template<
typename TypeTag>
1538 std::vector<double>& well_potentials,
1542 this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
1549 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1560 const auto&
ws = well_state.well(this->index_of_well_);
1561 if (this->isInjector())
1562 bhp = std::max(
ws.bhp, bhp);
1564 bhp = std::min(
ws.bhp, bhp);
1566 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1567 computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials,
deferred_logger);
1570 well_potentials = computeWellPotentialWithTHP(ebosSimulator,
deferred_logger, well_state);
1573 this->checkNegativeWellPotentials(well_potentials,
1574 this->param_.check_well_operability_,
1584 template<
typename TypeTag>
1599 template<
typename TypeTag>
1601 StandardWell<TypeTag>::
1602 updatePrimaryVariables(
const SummaryState& summary_state,
1603 const WellState& well_state,
1604 DeferredLogger& deferred_logger)
1606 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1608 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
1609 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1612 if constexpr (Base::has_polymermw) {
1613 this->primary_variables_.updatePolyMW(well_state);
1616 this->primary_variables_.checkFinite(deferred_logger);
1622 template<
typename TypeTag>
1624 StandardWell<TypeTag>::
1625 getRefDensity()
const
1627 return this->connections_.rho();
1633 template<
typename TypeTag>
1635 StandardWell<TypeTag>::
1636 updateWaterMobilityWithPolymer(
const Simulator& ebos_simulator,
1638 std::vector<EvalWell>& mob,
1639 DeferredLogger& deferred_logger)
const
1641 const int cell_idx = this->well_cells_[perf];
1642 const auto& int_quant = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1643 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1647 if (this->isInjector()) {
1649 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1650 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1651 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1654 if (PolymerModule::hasPlyshlog()) {
1657 if (this->isInjector() && this->wpolymer() == 0.) {
1662 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
1663 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1665 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1666 PerforationRates perf_rates;
1667 double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
1668 const double Tw = this->well_index_[perf] * trans_mult;
1669 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1670 perf_rates, deferred_logger);
1672 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1673 const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
1674 const auto& scaled_drainage_info =
1675 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1676 const double swcr = scaled_drainage_info.Swcr;
1677 const EvalWell poro = this->extendEval(int_quant.porosity());
1678 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1680 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1681 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1682 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1684 if (PolymerModule::hasShrate()) {
1687 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1689 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1690 int_quant.pvtRegionIndex(),
1693 mob[waterCompIdx] /= shear_factor;
1697 template<
typename TypeTag>
1699 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
1701 this->linSys_.extract(jacobian);
1705 template <
typename TypeTag>
1707 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1708 const BVector& weights,
1709 const int pressureVarIndex,
1710 const bool use_well_weights,
1711 const WellState& well_state)
const
1713 this->linSys_.extractCPRPressureMatrix(jacobian,
1724 template<
typename TypeTag>
1725 typename StandardWell<TypeTag>::EvalWell
1726 StandardWell<TypeTag>::
1727 pskinwater(
const double throughput,
1728 const EvalWell& water_velocity,
1729 DeferredLogger& deferred_logger)
const
1731 if constexpr (Base::has_polymermw) {
1732 const int water_table_id = this->polymerWaterTable_();
1733 if (water_table_id <= 0) {
1734 OPM_DEFLOG_THROW(std::runtime_error,
1735 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1738 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1739 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1741 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1742 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1745 OPM_DEFLOG_THROW(std::runtime_error,
1746 fmt::format(
"Polymermw is not activated, while injecting "
1747 "skin pressure is requested for well {}", name()),
1756 template<
typename TypeTag>
1757 typename StandardWell<TypeTag>::EvalWell
1758 StandardWell<TypeTag>::
1759 pskin(
const double throughput,
1760 const EvalWell& water_velocity,
1761 const EvalWell& poly_inj_conc,
1762 DeferredLogger& deferred_logger)
const
1764 if constexpr (Base::has_polymermw) {
1765 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
1766 const EvalWell water_velocity_abs = abs(water_velocity);
1767 if (poly_inj_conc == 0.) {
1768 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1770 const int polymer_table_id = this->polymerTable_();
1771 if (polymer_table_id <= 0) {
1772 OPM_DEFLOG_THROW(std::runtime_error,
1773 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1776 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1777 const double reference_concentration = skprpolytable.refConcentration;
1778 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1780 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1781 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1782 if (poly_inj_conc == reference_concentration) {
1783 return sign * pskin_poly;
1786 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1787 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1788 return sign * pskin;
1790 OPM_DEFLOG_THROW(std::runtime_error,
1791 fmt::format(
"Polymermw is not activated, while injecting "
1792 "skin pressure is requested for well {}", name()),
1801 template<
typename TypeTag>
1802 typename StandardWell<TypeTag>::EvalWell
1803 StandardWell<TypeTag>::
1804 wpolymermw(
const double throughput,
1805 const EvalWell& water_velocity,
1806 DeferredLogger& deferred_logger)
const
1808 if constexpr (Base::has_polymermw) {
1809 const int table_id = this->polymerInjTable_();
1810 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
1811 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1812 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
1813 if (this->wpolymer() == 0.) {
1814 return molecular_weight;
1816 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
1817 return molecular_weight;
1819 OPM_DEFLOG_THROW(std::runtime_error,
1820 fmt::format(
"Polymermw is not activated, while injecting "
1821 "polymer molecular weight is requested for well {}", name()),
1830 template<
typename TypeTag>
1832 StandardWell<TypeTag>::
1833 updateWaterThroughput(
const double dt, WellState &well_state)
const
1835 if constexpr (Base::has_polymermw) {
1836 if (this->isInjector()) {
1837 auto& ws = well_state.well(this->index_of_well_);
1838 auto& perf_water_throughput = ws.perf_data.water_throughput;
1839 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1840 const double perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
1842 if (perf_water_vel > 0.) {
1843 perf_water_throughput[perf] += perf_water_vel * dt;
1854 template<
typename TypeTag>
1856 StandardWell<TypeTag>::
1857 handleInjectivityRate(
const Simulator& ebosSimulator,
1859 std::vector<EvalWell>& cq_s)
const
1861 const int cell_idx = this->well_cells_[perf];
1862 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1863 const auto& fs = int_quants.fluidState();
1864 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
1865 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
1866 const int wat_vel_index = Bhp + 1 + perf;
1867 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1871 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
1877 template<
typename TypeTag>
1879 StandardWell<TypeTag>::
1880 handleInjectivityEquations(
const Simulator& ebosSimulator,
1881 const WellState& well_state,
1883 const EvalWell& water_flux_s,
1884 DeferredLogger& deferred_logger)
1886 const int cell_idx = this->well_cells_[perf];
1887 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1888 const auto& fs = int_quants.fluidState();
1889 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
1890 const EvalWell water_flux_r = water_flux_s / b_w;
1891 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
1892 const EvalWell water_velocity = water_flux_r / area;
1893 const int wat_vel_index = Bhp + 1 + perf;
1896 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
1898 const auto& ws = well_state.well(this->index_of_well_);
1899 const auto& perf_data = ws.perf_data;
1900 const auto& perf_water_throughput = perf_data.water_throughput;
1901 const double throughput = perf_water_throughput[perf];
1902 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
1904 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1905 poly_conc.setValue(this->wpolymer());
1908 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
1909 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
1911 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
1912 assembleInjectivityEq(eq_pskin,
1917 this->primary_variables_.numWellEq(),
1925 template<
typename TypeTag>
1927 StandardWell<TypeTag>::
1928 checkConvergenceExtraEqs(
const std::vector<double>& res,
1929 ConvergenceReport& report)
const
1934 if constexpr (Base::has_polymermw) {
1935 WellConvergence(*this).
1936 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
1944 template<
typename TypeTag>
1946 StandardWell<TypeTag>::
1947 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
1948 const IntensiveQuantities& int_quants,
1949 const WellState& well_state,
1951 std::vector<RateVector>& connectionRates,
1952 DeferredLogger& deferred_logger)
const
1955 EvalWell cq_s_polymw = cq_s_poly;
1956 if (this->isInjector()) {
1957 const int wat_vel_index = Bhp + 1 + perf;
1958 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
1959 if (water_velocity > 0.) {
1960 const auto& ws = well_state.well(this->index_of_well_);
1961 const auto& perf_water_throughput = ws.perf_data.water_throughput;
1962 const double throughput = perf_water_throughput[perf];
1963 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
1964 cq_s_polymw *= molecular_weight;
1970 }
else if (this->isProducer()) {
1971 if (cq_s_polymw < 0.) {
1972 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
1979 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
1987 template<
typename TypeTag>
1988 std::optional<double>
1989 StandardWell<TypeTag>::
1990 computeBhpAtThpLimitProd(
const WellState& well_state,
1991 const Simulator& ebos_simulator,
1992 const SummaryState& summary_state,
1993 DeferredLogger& deferred_logger)
const
1995 return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
1997 this->getALQ(well_state),
2001 template<
typename TypeTag>
2002 std::optional<double>
2003 StandardWell<TypeTag>::
2004 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
2005 const SummaryState& summary_state,
2006 const double alq_value,
2007 DeferredLogger& deferred_logger)
const
2010 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2016 std::vector<double> rates(3);
2017 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2018 this->adaptRatesForVFP(rates);
2022 double max_pressure = 0.0;
2023 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2024 const int cell_idx = this->well_cells_[perf];
2025 const auto& int_quants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
2026 const auto& fs = int_quants.fluidState();
2027 double pressure_cell = this->getPerfCellPressure(fs).value();
2028 max_pressure = std::max(max_pressure, pressure_cell);
2030 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2033 this->connections_.rho(),
2035 this->getTHPConstraint(summary_state),
2039 auto v = frates(*bhpAtLimit);
2040 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2046 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2050 std::vector<double> rates(3);
2051 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
2052 this->adaptRatesForVFP(rates);
2056 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2059 this->connections_.rho(),
2061 this->getTHPConstraint(summary_state),
2067 auto v = frates(*bhpAtLimit);
2068 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2074 return std::nullopt;
2079 template<
typename TypeTag>
2080 std::optional<double>
2081 StandardWell<TypeTag>::
2082 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
2083 const SummaryState& summary_state,
2084 DeferredLogger& deferred_logger)
const
2087 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2093 std::vector<double> rates(3);
2094 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2098 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2100 this->connections_.rho(),
2111 template<
typename TypeTag>
2113 StandardWell<TypeTag>::
2114 iterateWellEqWithControl(
const Simulator& ebosSimulator,
2116 const Well::InjectionControls& inj_controls,
2117 const Well::ProductionControls& prod_controls,
2118 WellState& well_state,
2119 const GroupState& group_state,
2120 DeferredLogger& deferred_logger)
2122 const int max_iter = this->param_.max_inner_iter_wells_;
2125 bool relax_convergence =
false;
2126 this->regularize_ =
false;
2127 const auto& summary_state = ebosSimulator.vanguard().summaryState();
2129 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2131 if (it > this->param_.strict_inner_iter_wells_) {
2132 relax_convergence =
true;
2133 this->regularize_ =
true;
2136 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2138 converged = report.converged();
2144 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2151 initPrimaryVariablesEvaluation();
2152 }
while (it < max_iter);
2158 template<
typename TypeTag>
2160 StandardWell<TypeTag>::
2161 iterateWellEqWithSwitching(
const Simulator& ebosSimulator,
2163 const Well::InjectionControls& inj_controls,
2164 const Well::ProductionControls& prod_controls,
2165 WellState& well_state,
2166 const GroupState& group_state,
2167 DeferredLogger& deferred_logger)
2169 const int max_iter = this->param_.max_inner_iter_wells_;
2173 bool relax_convergence =
false;
2174 this->regularize_ =
false;
2175 const auto& summary_state = ebosSimulator.vanguard().summaryState();
2178 constexpr int min_its_after_switch = 2;
2179 int its_since_last_switch = min_its_after_switch;
2180 int switch_count= 0;
2181 const auto well_status = this->wellStatus_;
2182 const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
2183 bool changed =
false;
2184 bool final_check =
false;
2186 its_since_last_switch++;
2187 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2188 const double wqTotal = this->primary_variables_.eval(WQTotal).value();
2189 changed = this->updateWellControlAndStatusLocalIteration(ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger);
2191 its_since_last_switch = 0;
2194 if (!changed && final_check) {
2197 final_check =
false;
2201 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2203 if (it > this->param_.strict_inner_iter_wells_) {
2204 relax_convergence =
true;
2205 this->regularize_ =
true;
2208 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2210 converged = report.converged();
2214 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2216 its_since_last_switch = min_its_after_switch;
2223 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2224 initPrimaryVariablesEvaluation();
2225 }
while (it < max_iter);
2228 if (allow_switching){
2230 const bool is_stopped = this->wellIsStopped();
2231 if (this->wellHasTHPConstraints(summary_state)){
2232 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2233 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2235 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2244 this->wellStatus_ = well_status;
2247 const std::string message = fmt::format(
" Well {} did not converged in {} inner iterations ("
2248 "{} control/status switches).", this->name(), it, switch_count);
2249 deferred_logger.debug(message);
2255 template<
typename TypeTag>
2262 std::vector<double>
well_q_s(this->num_components_, 0.);
2263 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2264 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2265 for (
int perf = 0;
perf < this->number_of_perforations_; ++
perf) {
2268 std::vector<Scalar>
mob(this->num_components_, 0.);
2270 std::vector<Scalar>
cq_s(this->num_components_, 0.);
2276 for (
int comp = 0;
comp < this->num_components_; ++
comp) {
2280 const auto& comm = this->parallel_well_info_.communication();
2281 if (comm.size() > 1)
2290 template <
typename TypeTag>
2295 const int num_pri_vars = this->primary_variables_.numWellEq();
2298 retval[
ii] = this->primary_variables_.value(
ii);
2307 template <
typename TypeTag>
2309 StandardWell<TypeTag>::
2310 setPrimaryVars(std::vector<double>::const_iterator it)
2312 const int num_pri_vars = this->primary_variables_.numWellEq();
2313 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2314 this->primary_variables_.setValue(ii, it[ii]);
2316 return num_pri_vars;
2320 template <
typename TypeTag>
2321 typename StandardWell<TypeTag>::Eval
2322 StandardWell<TypeTag>::
2323 connectionRateEnergy(
const double maxOilSaturation,
2324 const std::vector<EvalWell>& cq_s,
2325 const IntensiveQuantities& intQuants,
2326 DeferredLogger& deferred_logger)
const
2328 auto fs = intQuants.fluidState();
2330 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2331 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2336 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2337 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2338 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2339 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2340 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2343 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2344 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2349 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2351 deferred_logger.debug(
2352 fmt::format(
"Problematic d value {} obtained for well {}"
2353 " during calculateSinglePerf with rs {}"
2354 ", rv {}. Continue as if no dissolution (rs = 0) and"
2355 " vaporization (rv = 0) for this connection.",
2356 d, this->name(), fs.Rs(), fs.Rv()));
2357 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2359 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2360 cq_r_thermal = (cq_s[gasCompIdx] -
2361 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2362 (d * this->extendEval(fs.invB(phaseIdx)) );
2363 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2365 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2367 (d * this->extendEval(fs.invB(phaseIdx)) );
2373 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
2375 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2376 fs.setTemperature(this->well_ecl_.temperature());
2377 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2378 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2379 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2380 paramCache.setRegionIndex(pvtRegionIdx);
2381 paramCache.setMaxOilSat(maxOilSaturation);
2382 paramCache.updatePhase(fs, phaseIdx);
2384 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2385 fs.setDensity(phaseIdx, rho);
2386 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2387 fs.setEnthalpy(phaseIdx, h);
2388 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2389 result += getValue(cq_r_thermal);
2392 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2393 result += Base::restrictEval(cq_r_thermal);
2401 template <
typename TypeTag>
2402 template<
class Value>
2404 StandardWell<TypeTag>::
2405 gasOilPerfRateInj(
const std::vector<Value>& cq_s,
2406 PerforationRates& perf_rates,
2409 const Value& pressure,
2411 DeferredLogger& deferred_logger)
const
2413 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2414 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2423 const double d = 1.0 - getValue(rv) * getValue(rs);
2426 deferred_logger.debug(dValueError(d, this->name(),
2427 "gasOilPerfRateInj",
2432 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2435 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
2438 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2443 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2449 template <
typename TypeTag>
2450 template<
class Value>
2452 StandardWell<TypeTag>::
2453 gasOilPerfRateProd(std::vector<Value>& cq_s,
2454 PerforationRates& perf_rates,
2457 const Value& rvw)
const
2459 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2460 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2461 const Value cq_sOil = cq_s[oilCompIdx];
2462 const Value cq_sGas = cq_s[gasCompIdx];
2463 const Value dis_gas = rs * cq_sOil;
2464 const Value vap_oil = rv * cq_sGas;
2466 cq_s[gasCompIdx] += dis_gas;
2467 cq_s[oilCompIdx] += vap_oil;
2470 if (this->isProducer()) {
2471 perf_rates.dis_gas = getValue(dis_gas);
2472 perf_rates.vap_oil = getValue(vap_oil);
2475 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2476 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2477 const Value vap_wat = rvw * cq_sGas;
2478 cq_s[waterCompIdx] += vap_wat;
2479 if (this->isProducer())
2480 perf_rates.vap_wat = getValue(vap_wat);
2485 template <
typename TypeTag>
2486 template<
class Value>
2488 StandardWell<TypeTag>::
2489 gasWaterPerfRateProd(std::vector<Value>& cq_s,
2490 PerforationRates& perf_rates,
2492 const Value& rsw)
const
2494 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2495 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2496 const Value cq_sWat = cq_s[waterCompIdx];
2497 const Value cq_sGas = cq_s[gasCompIdx];
2498 const Value vap_wat = rvw * cq_sGas;
2499 const Value dis_gas_wat = rsw * cq_sWat;
2500 cq_s[waterCompIdx] += vap_wat;
2501 cq_s[gasCompIdx] += dis_gas_wat;
2502 if (this->isProducer()) {
2503 perf_rates.vap_wat = getValue(vap_wat);
2504 perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
2509 template <
typename TypeTag>
2510 template<
class Value>
2512 StandardWell<TypeTag>::
2513 gasWaterPerfRateInj(
const std::vector<Value>& cq_s,
2514 PerforationRates& perf_rates,
2517 const Value& pressure,
2518 DeferredLogger& deferred_logger)
const
2521 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2522 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2524 const double dw = 1.0 - getValue(rvw) * getValue(rsw);
2527 deferred_logger.debug(dValueError(dw, this->name(),
2528 "gasWaterPerfRateInj",
2529 rsw, rvw, pressure));
2533 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
2536 perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
2541 template <
typename TypeTag>
2542 template<
class Value>
2544 StandardWell<TypeTag>::
2545 disOilVapWatVolumeRatio(Value& volumeRatio,
2548 const Value& pressure,
2549 const std::vector<Value>& cmix_s,
2550 const std::vector<Value>& b_perfcells_dense,
2551 DeferredLogger& deferred_logger)
const
2553 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2554 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2556 const Value d = 1.0 - rvw * rsw;
2559 deferred_logger.debug(dValueError(d, this->name(),
2560 "disOilVapWatVolumeRatio",
2561 rsw, rvw, pressure));
2563 const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
2564 : cmix_s[waterCompIdx];
2565 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
2567 const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
2568 : cmix_s[waterCompIdx];
2569 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
2573 template <
typename TypeTag>
2574 template<
class Value>
2576 StandardWell<TypeTag>::
2577 gasOilVolumeRatio(Value& volumeRatio,
2580 const Value& pressure,
2581 const std::vector<Value>& cmix_s,
2582 const std::vector<Value>& b_perfcells_dense,
2583 DeferredLogger& deferred_logger)
const
2585 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2586 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2588 const Value d = 1.0 - rv * rs;
2591 deferred_logger.debug(dValueError(d, this->name(),
2592 "gasOilVolumeRatio",
2595 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
2596 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
2598 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
2599 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
Definition AquiferInterface.hpp:35
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
Definition StandardWell.hpp:60
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition StandardWell_impl.hpp:1295
Class for computing BHP limits.
Definition WellBhpThpCalculator.hpp:42
double mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
Definition WellBhpThpCalculator.cpp:89
Collect per-connection static information to enable calculating connection-level or well-level produc...
Definition WellProdIndexCalculator.hpp:36
double connectionProdIndStandard(const std::size_t connIdx, const double connMobility) const
Compute connection-level steady-state productivity index value using dynamic phase mobility.
Definition WellProdIndexCalculator.cpp:110
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Definition PerforationData.hpp:39