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>&>(*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,
122 const std::vector<Scalar>& Tw,
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,
227 const std::vector<Scalar>& Tw,
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[componentIdx] * (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 Value volumeRatio = bhp * 0.0;
276 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
277 disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
278 cmix_s, b_perfcells_dense, deferred_logger);
280 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
281 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
282 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
286 if constexpr (Indices::enableSolvent) {
287 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
290 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
291 gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
292 cmix_s, b_perfcells_dense, deferred_logger);
295 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
296 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
297 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
299 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
300 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
301 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
306 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
307 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
308 Value cqt_is = cqt_i / volumeRatio;
309 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
313 if (this->isProducer()) {
314 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
315 gasOilPerfRateInj(cq_s, perf_rates,
316 rv, rs, pressure, rvw, deferred_logger);
318 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
320 gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
321 pressure, deferred_logger);
328 template<
typename TypeTag>
330 StandardWell<TypeTag>::
331 assembleWellEqWithoutIteration(
const Simulator& simulator,
333 const Well::InjectionControls& inj_controls,
334 const Well::ProductionControls& prod_controls,
335 WellState& well_state,
336 const GroupState& group_state,
337 DeferredLogger& deferred_logger)
341 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
344 this->linSys_.clear();
346 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
347 prod_controls, well_state,
348 group_state, deferred_logger);
354 template<
typename TypeTag>
356 StandardWell<TypeTag>::
357 assembleWellEqWithoutIterationImpl(
const Simulator& simulator,
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);
373 const int np = this->number_of_phases_;
375 std::vector<RateVector> connectionRates = this->connectionRates_;
377 auto& perf_data = ws.perf_data;
378 auto& perf_rates = perf_data.phase_rates;
379 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
381 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
382 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
383 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
384 calculateSinglePerf(simulator, perf, well_state, connectionRates,
385 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
388 if constexpr (has_polymer && Base::has_polymermw) {
389 if (this->isInjector()) {
390 handleInjectivityEquations(simulator, well_state, perf,
391 water_flux_s, deferred_logger);
394 const int cell_idx = this->well_cells_[perf];
395 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
397 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
399 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
401 StandardWellAssemble<FluidSystem,Indices>(*this).
402 assemblePerforationEq(cq_s_effective,
405 this->primary_variables_.numWellEq(),
409 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
410 auto& perf_rate_solvent = perf_data.solvent_rates;
411 perf_rate_solvent[perf] = cq_s[componentIdx].value();
413 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
417 if constexpr (has_zFraction) {
418 StandardWellAssemble<FluidSystem,Indices>(*this).
419 assembleZFracEq(cq_s_zfrac_effective,
421 this->primary_variables_.numWellEq(),
426 this->connectionRates_ = connectionRates;
431 const auto& comm = this->parallel_well_info_.communication();
432 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
436 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
439 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
442 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
443 if (FluidSystem::numActivePhases() > 1) {
445 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
446 this->F0_[componentIdx]) * volume / dt;
448 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
449 StandardWellAssemble<FluidSystem,Indices>(*this).
450 assembleSourceEq(resWell_loc,
452 this->primary_variables_.numWellEq(),
456 const auto& summaryState = simulator.vanguard().summaryState();
457 const Schedule& schedule = simulator.vanguard().schedule();
458 StandardWellAssemble<FluidSystem,Indices>(*this).
459 assembleControlEq(well_state, group_state,
460 schedule, summaryState,
461 inj_controls, prod_controls,
462 this->primary_variables_,
463 this->connections_.rho(),
470 this->linSys_.invert();
472 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
479 template<
typename TypeTag>
481 StandardWell<TypeTag>::
482 calculateSinglePerf(
const Simulator& simulator,
484 WellState& well_state,
485 std::vector<RateVector>& connectionRates,
486 std::vector<EvalWell>& cq_s,
487 EvalWell& water_flux_s,
488 EvalWell& cq_s_zfrac_effective,
489 DeferredLogger& deferred_logger)
const
491 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
492 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
493 const int cell_idx = this->well_cells_[perf];
494 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
495 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
496 getMobility(simulator, perf, mob, deferred_logger);
498 PerforationRates perf_rates;
499 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
500 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
501 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
502 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
503 cq_s, perf_rates, deferred_logger);
505 auto& ws = well_state.well(this->index_of_well_);
506 auto& perf_data = ws.perf_data;
507 if constexpr (has_polymer && Base::has_polymermw) {
508 if (this->isInjector()) {
511 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
512 water_flux_s = cq_s[water_comp_idx];
515 handleInjectivityRate(simulator, perf, cq_s);
520 if (this->isProducer()) {
521 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
522 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
523 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
524 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
525 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.dis_gas;
526 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.dis_gas_in_water;
527 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.vap_oil;
528 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.vap_wat;
531 if constexpr (has_energy) {
532 connectionRates[perf][Indices::contiEnergyEqIdx] =
533 connectionRateEnergy(simulator.problem().maxOilSaturation(cell_idx),
534 cq_s, intQuants, deferred_logger);
537 if constexpr (has_polymer) {
538 std::variant<Scalar,EvalWell> polymerConcentration;
539 if (this->isInjector()) {
540 polymerConcentration = this->wpolymer();
542 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
543 intQuants.polymerViscosityCorrection());
546 [[maybe_unused]] EvalWell cq_s_poly;
547 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
549 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
550 cq_s, polymerConcentration);
552 if constexpr (Base::has_polymermw) {
553 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
554 perf, connectionRates, deferred_logger);
558 if constexpr (has_foam) {
559 std::variant<Scalar,EvalWell> foamConcentration;
560 if (this->isInjector()) {
561 foamConcentration = this->wfoam();
563 foamConcentration = this->extendEval(intQuants.foamConcentration());
565 connectionRates[perf][Indices::contiFoamEqIdx] =
566 this->connections_.connectionRateFoam(cq_s, foamConcentration,
567 FoamModule::transportPhase(),
571 if constexpr (has_zFraction) {
572 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
573 if (this->isInjector()) {
574 solventConcentration = this->wsolvent();
576 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
577 this->extendEval(intQuants.yVolume())};
579 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
580 cq_s_zfrac_effective) =
581 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
582 perf_rates.dis_gas, cq_s,
583 solventConcentration);
586 if constexpr (has_brine) {
587 std::variant<Scalar,EvalWell> saltConcentration;
588 if (this->isInjector()) {
589 saltConcentration = this->wsalt();
591 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
594 connectionRates[perf][Indices::contiBrineEqIdx] =
595 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
596 perf_rates.vap_wat, cq_s,
600 if constexpr (has_micp) {
601 std::variant<Scalar,EvalWell> microbialConcentration;
602 std::variant<Scalar,EvalWell> oxygenConcentration;
603 std::variant<Scalar,EvalWell> ureaConcentration;
604 if (this->isInjector()) {
605 microbialConcentration = this->wmicrobes();
606 oxygenConcentration = this->woxygen();
607 ureaConcentration = this->wurea();
609 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
610 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
611 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
613 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
614 connectionRates[perf][Indices::contiOxygenEqIdx],
615 connectionRates[perf][Indices::contiUreaEqIdx]) =
616 this->connections_.connectionRatesMICP(cq_s,
617 microbialConcentration,
623 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
628 template<
typename TypeTag>
629 template<
class Value>
631 StandardWell<TypeTag>::
632 getMobility(
const Simulator& simulator,
634 std::vector<Value>& mob,
635 DeferredLogger& deferred_logger)
const
637 auto obtain = [
this](
const Eval& value)
639 if constexpr (std::is_same_v<Value, Scalar>) {
640 static_cast<void>(
this);
641 return getValue(value);
643 return this->extendEval(value);
646 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
647 obtain, deferred_logger);
650 if constexpr (has_polymer) {
651 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
652 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
657 if constexpr (!Base::has_polymermw) {
658 if constexpr (std::is_same_v<Value, Scalar>) {
659 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
660 for (std::size_t i = 0; i < mob.size(); ++i) {
661 mob_eval[i].setValue(mob[i]);
663 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
664 for (std::size_t i = 0; i < mob.size(); ++i) {
665 mob[i] = getValue(mob_eval[i]);
668 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
674 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
675 const double bhp = this->primary_variables_.value(Bhp);
676 const double perf_press = bhp + this->connections_.pressure_diff(perf);
677 const double multiplier = this->getInjMult(perf, bhp, perf_press);
678 for (std::size_t i = 0; i < mob.size(); ++i) {
679 mob[i] *= multiplier;
685 template<
typename TypeTag>
687 StandardWell<TypeTag>::
688 updateWellState(
const SummaryState& summary_state,
689 const BVectorWell& dwells,
690 WellState& well_state,
691 DeferredLogger& deferred_logger)
693 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
695 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
696 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
698 updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
699 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
706 template<
typename TypeTag>
708 StandardWell<TypeTag>::
709 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
710 const bool stop_or_zero_rate_target,
711 DeferredLogger& deferred_logger)
713 const double dFLimit = this->param_.dwell_fraction_max_;
714 const double dBHPLimit = this->param_.dbhp_max_rel_;
715 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
718 if constexpr (Base::has_polymermw) {
719 this->primary_variables_.updateNewtonPolyMW(dwells);
722 this->primary_variables_.checkFinite(deferred_logger);
729 template<
typename TypeTag>
731 StandardWell<TypeTag>::
732 updateWellStateFromPrimaryVariables(
const bool stop_or_zero_rate_target,
733 WellState& well_state,
734 const SummaryState& summary_state,
735 DeferredLogger& deferred_logger)
const
737 this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
740 if constexpr (Base::has_polymermw) {
741 this->primary_variables_.copyToWellStatePolyMW(well_state);
749 template<
typename TypeTag>
751 StandardWell<TypeTag>::
752 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
757 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
758 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
760 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
761 std::vector<Scalar> mob(this->num_components_, 0.0);
762 getMobility(simulator, perf, mob, deferred_logger);
764 const int cell_idx = this->well_cells_[perf];
765 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
766 const auto& fs = int_quantities.fluidState();
768 double p_r = this->getPerfCellPressure(fs).value();
771 std::vector<double> b_perf(this->num_components_);
772 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
773 if (!FluidSystem::phaseIsActive(phase)) {
776 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
777 b_perf[comp_idx] = fs.invB(phase).value();
779 if constexpr (has_solvent) {
780 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
784 const double h_perf = this->connections_.pressure_diff(perf);
785 const double pressure_diff = p_r - h_perf;
790 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
791 deferred_logger.debug(
"CROSSFLOW_IPR",
792 "cross flow found when updateIPR for well " + name()
793 +
" . The connection is ignored in IPR calculations");
799 double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
800 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
801 const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
802 std::vector<double> ipr_a_perf(this->ipr_a_.size());
803 std::vector<double> ipr_b_perf(this->ipr_b_.size());
804 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
805 const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
806 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
807 ipr_b_perf[comp_idx] += tw_mob;
811 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
812 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
813 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
814 const double rs = (fs.Rs()).value();
815 const double rv = (fs.Rv()).value();
817 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
818 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
820 ipr_a_perf[gas_comp_idx] += dis_gas_a;
821 ipr_a_perf[oil_comp_idx] += vap_oil_a;
823 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
824 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
826 ipr_b_perf[gas_comp_idx] += dis_gas_b;
827 ipr_b_perf[oil_comp_idx] += vap_oil_b;
830 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
831 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
832 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
835 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
836 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
839 template<
typename TypeTag>
841 StandardWell<TypeTag>::
842 updateIPRImplicit(
const Simulator& simulator,
843 WellState& well_state,
844 DeferredLogger& deferred_logger)
853 auto rates = well_state.well(this->index_of_well_).surface_rates;
855 for (std::size_t p = 0; p < rates.size(); ++p) {
856 zero_rates &= rates[p] == 0.0;
858 auto& ws = well_state.well(this->index_of_well_);
860 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
861 deferred_logger.debug(msg);
873 const auto& group_state = simulator.problem().wellModel().groupState();
875 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
876 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
878 auto inj_controls = Well::InjectionControls(0);
879 auto prod_controls = Well::ProductionControls(0);
880 prod_controls.addControl(Well::ProducerCMode::BHP);
881 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
884 const auto cmode = ws.production_cmode;
885 ws.production_cmode = Well::ProducerCMode::BHP;
886 const double dt = simulator.timeStepSize();
887 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
889 const size_t nEq = this->primary_variables_.numWellEq();
893 for (
size_t i=0; i < nEq; ++i){
898 BVectorWell x_well(1);
899 x_well[0].resize(nEq);
900 this->linSys_.solve(rhs, x_well);
902 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
903 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
904 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
905 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
907 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
909 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
912 ws.production_cmode = cmode;
915 template<
typename TypeTag>
917 StandardWell<TypeTag>::
918 checkOperabilityUnderBHPLimit(
const WellState& well_state,
919 const Simulator& simulator,
920 DeferredLogger& deferred_logger)
922 const auto& summaryState = simulator.vanguard().summaryState();
923 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
926 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
927 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
930 double total_ipr_mass_rate = 0.0;
931 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
933 if (!FluidSystem::phaseIsActive(phaseIdx)) {
937 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
938 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
940 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
941 total_ipr_mass_rate += ipr_rate * rho;
943 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
944 this->operability_status_.operable_under_only_bhp_limit =
false;
948 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
952 std::vector<double> well_rates_bhp_limit;
953 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
955 this->adaptRatesForVFP(well_rates_bhp_limit);
956 const double thp_limit = this->getTHPConstraint(summaryState);
957 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
959 this->connections_.rho(),
960 this->getALQ(well_state),
963 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
964 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
975 this->operability_status_.operable_under_only_bhp_limit =
true;
976 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
984 template<
typename TypeTag>
986 StandardWell<TypeTag>::
987 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
988 const WellState& well_state,
989 DeferredLogger& deferred_logger)
991 const auto& summaryState = simulator.vanguard().summaryState();
992 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
993 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
996 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
998 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
999 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1000 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1002 const double thp_limit = this->getTHPConstraint(summaryState);
1003 if (this->isProducer() && *obtain_bhp < thp_limit) {
1004 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1005 +
" bars is SMALLER than thp limit "
1006 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1007 +
" bars as a producer for well " + name();
1008 deferred_logger.debug(msg);
1010 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1011 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1012 +
" bars is LARGER than thp limit "
1013 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1014 +
" bars as a injector for well " + name();
1015 deferred_logger.debug(msg);
1018 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1019 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1020 if (!this->wellIsStopped()) {
1021 const double thp_limit = this->getTHPConstraint(summaryState);
1022 deferred_logger.debug(
" could not find bhp value at thp limit "
1023 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1024 +
" bar for well " + name() +
", the well might need to be closed ");
1033 template<
typename TypeTag>
1035 StandardWell<TypeTag>::
1036 allDrawDownWrongDirection(
const Simulator& simulator)
const
1038 bool all_drawdown_wrong_direction =
true;
1040 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1041 const int cell_idx = this->well_cells_[perf];
1042 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1043 const auto& fs = intQuants.fluidState();
1045 const double pressure = this->getPerfCellPressure(fs).value();
1046 const double bhp = this->primary_variables_.eval(Bhp).value();
1049 const double well_pressure = bhp + this->connections_.pressure_diff(perf);
1050 const double drawdown = pressure - well_pressure;
1055 if ( (drawdown < 0. && this->isInjector()) ||
1056 (drawdown > 0. && this->isProducer()) ) {
1057 all_drawdown_wrong_direction =
false;
1062 const auto& comm = this->parallel_well_info_.communication();
1063 if (comm.size() > 1)
1065 all_drawdown_wrong_direction =
1066 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1069 return all_drawdown_wrong_direction;
1075 template<
typename TypeTag>
1077 StandardWell<TypeTag>::
1078 canProduceInjectWithCurrentBhp(
const Simulator& simulator,
1079 const WellState& well_state,
1080 DeferredLogger& deferred_logger)
1082 const double bhp = well_state.well(this->index_of_well_).bhp;
1083 std::vector<double> well_rates;
1084 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1086 const double sign = (this->isProducer()) ? -1. : 1.;
1087 const double threshold = sign * std::numeric_limits<double>::min();
1089 bool can_produce_inject =
false;
1090 for (
const auto value : well_rates) {
1091 if (this->isProducer() && value < threshold) {
1092 can_produce_inject =
true;
1094 }
else if (this->isInjector() && value > threshold) {
1095 can_produce_inject =
true;
1100 if (!can_produce_inject) {
1101 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1104 return can_produce_inject;
1111 template<
typename TypeTag>
1113 StandardWell<TypeTag>::
1114 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1116 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1122 template<
typename TypeTag>
1124 StandardWell<TypeTag>::
1125 computePropertiesForWellConnectionPressures(
const Simulator& simulator,
1126 const WellState& well_state,
1127 WellConnectionProps& props)
const
1129 std::function<Scalar(
int,
int)> getTemperature =
1130 [&simulator](
int cell_idx,
int phase_idx)
1132 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().temperature(phase_idx).value();
1134 std::function<Scalar(
int)> getSaltConcentration =
1135 [&simulator](
int cell_idx)
1137 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().saltConcentration().value();
1139 std::function<int(
int)> getPvtRegionIdx =
1140 [&simulator](
int cell_idx)
1142 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().pvtRegionIndex();
1144 std::function<Scalar(
int)> getInvFac =
1145 [&simulator](
int cell_idx)
1147 return simulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1149 std::function<Scalar(
int)> getSolventDensity =
1150 [&simulator](
int cell_idx)
1152 return simulator.model().intensiveQuantities(cell_idx, 0).solventRefDensity();
1155 this->connections_.computePropertiesForPressures(well_state,
1157 getSaltConcentration,
1168 template<
typename TypeTag>
1173 const std::vector<double>& B_avg,
1175 const bool relax_tolerance)
const
1179 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1181 double tol_wells = this->param_.tolerance_wells_;
1183 constexpr double stopped_factor = 1.e-4;
1185 constexpr double dynamic_thp_factor = 1.e-1;
1186 if (this->stopppedOrZeroRateTarget(summary_state, well_state)) {
1187 tol_wells = tol_wells*stopped_factor;
1188 }
else if (this->getDynamicThpLimit()) {
1189 tol_wells = tol_wells*dynamic_thp_factor;
1192 std::vector<double> res;
1195 this->param_.max_residual_allowed_,
1197 this->param_.relaxed_tolerance_flow_well_,
1199 this->wellIsStopped(),
1203 checkConvergenceExtraEqs(res, report);
1212 template<
typename TypeTag>
1220 auto fluidState = [&simulator,
this](
const int perf)
1222 const auto cell_idx = this->well_cells_[perf];
1223 return simulator.model()
1224 .intensiveQuantities(cell_idx, 0).fluidState();
1227 const int np = this->number_of_phases_;
1228 auto setToZero = [np](
double* x) ->
void
1230 std::fill_n(x, np, 0.0);
1233 auto addVector = [np](
const double* src,
double* dest) ->
void
1235 std::transform(src, src + np, dest, dest, std::plus<>{});
1238 auto& ws = well_state.well(this->index_of_well_);
1239 auto& perf_data = ws.perf_data;
1240 auto* wellPI = ws.productivity_index.data();
1241 auto* connPI = perf_data.prod_index.data();
1245 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1246 auto subsetPerfID = 0;
1248 for (
const auto& perf : *this->perf_data_) {
1249 auto allPerfID = perf.ecl_index;
1251 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1256 std::vector<Scalar> mob(this->num_components_, 0.0);
1257 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1259 const auto& fs = fluidState(subsetPerfID);
1262 if (this->isInjector()) {
1263 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1264 mob, connPI, deferred_logger);
1267 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1270 addVector(connPI, wellPI);
1277 const auto& comm = this->parallel_well_info_.communication();
1278 if (comm.size() > 1) {
1279 comm.sum(wellPI, np);
1282 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1283 "Internal logic error in processing connections for PI/II");
1288 template<
typename TypeTag>
1290 StandardWell<TypeTag>::
1291 computeWellConnectionDensitesPressures(
const Simulator& simulator,
1292 const WellState& well_state,
1293 const WellConnectionProps& props,
1294 DeferredLogger& deferred_logger)
1296 std::function<Scalar(
int,
int)> invB =
1297 [&simulator](
int cell_idx,
int phase_idx)
1299 return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().invB(phase_idx).value();
1301 std::function<Scalar(
int,
int)> mobility =
1302 [&simulator](
int cell_idx,
int phase_idx)
1304 return simulator.model().intensiveQuantities(cell_idx, 0).mobility(phase_idx).value();
1306 std::function<Scalar(
int)> invFac =
1307 [&simulator](
int cell_idx)
1309 return simulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1311 std::function<Scalar(
int)> solventMobility =
1312 [&simulator](
int cell_idx)
1314 return simulator.model().intensiveQuantities(cell_idx, 0).solventMobility().value();
1317 this->connections_.computeProperties(well_state,
1330 template<
typename TypeTag>
1332 StandardWell<TypeTag>::
1333 computeWellConnectionPressures(
const Simulator& simulator,
1334 const WellState& well_state,
1335 DeferredLogger& deferred_logger)
1340 WellConnectionProps props;
1341 computePropertiesForWellConnectionPressures(simulator, well_state, props);
1342 computeWellConnectionDensitesPressures(simulator, well_state,
1343 props, deferred_logger);
1350 template<
typename TypeTag>
1352 StandardWell<TypeTag>::
1353 solveEqAndUpdateWellState(
const SummaryState& summary_state,
1354 WellState& well_state,
1355 DeferredLogger& deferred_logger)
1357 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1361 BVectorWell dx_well(1);
1362 dx_well[0].resize(this->primary_variables_.numWellEq());
1363 this->linSys_.solve( dx_well);
1365 updateWellState(summary_state, dx_well, well_state, deferred_logger);
1372 template<
typename TypeTag>
1374 StandardWell<TypeTag>::
1375 calculateExplicitQuantities(
const Simulator& simulator,
1376 const WellState& well_state,
1377 DeferredLogger& deferred_logger)
1379 const auto& summary_state = simulator.vanguard().summaryState();
1380 updatePrimaryVariables(summary_state, well_state, deferred_logger);
1381 initPrimaryVariablesEvaluation();
1382 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1383 this->computeAccumWell();
1388 template<
typename TypeTag>
1391 apply(
const BVector& x, BVector& Ax)
const
1393 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1395 if (this->param_.matrix_add_well_contributions_)
1401 this->linSys_.
apply(x, Ax);
1407 template<
typename TypeTag>
1410 apply(BVector& r)
const
1412 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1414 this->linSys_.
apply(r);
1420 template<
typename TypeTag>
1428 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1431 xw[0].resize(this->primary_variables_.numWellEq());
1433 this->linSys_.recoverSolutionWell(x, xw);
1434 updateWellState(summary_state, xw, well_state, deferred_logger);
1440 template<
typename TypeTag>
1445 std::vector<double>& well_flux,
1449 const int np = this->number_of_phases_;
1450 well_flux.resize(np, 0.0);
1452 const bool allow_cf = this->getAllowCrossFlow();
1454 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1455 const int cell_idx = this->well_cells_[perf];
1456 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1458 std::vector<Scalar> mob(this->num_components_, 0.);
1459 getMobility(simulator, perf, mob, deferred_logger);
1460 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
1461 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1462 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1464 std::vector<Scalar> cq_s(this->num_components_, 0.);
1466 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1467 cq_s, perf_rates, deferred_logger);
1469 for(
int p = 0; p < np; ++p) {
1470 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1474 if constexpr (has_solvent) {
1475 const auto& pu = this->phaseUsage();
1476 assert(pu.phase_used[Gas]);
1477 const int gas_pos = pu.phase_pos[Gas];
1478 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1481 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1486 template<
typename TypeTag>
1488 StandardWell<TypeTag>::
1489 computeWellRatesWithBhpIterations(
const Simulator& simulator,
1491 std::vector<double>& well_flux,
1492 DeferredLogger& deferred_logger)
const
1496 StandardWell<TypeTag> well_copy(*
this);
1501 WellState well_state_copy = simulator.problem().wellModel().wellState();
1502 const auto& group_state = simulator.problem().wellModel().groupState();
1505 const auto& summary_state = simulator.vanguard().summaryState();
1506 auto inj_controls = well_copy.well_ecl_.isInjector()
1507 ? well_copy.well_ecl_.injectionControls(summary_state)
1508 : Well::InjectionControls(0);
1509 auto prod_controls = well_copy.well_ecl_.isProducer()
1510 ? well_copy.well_ecl_.productionControls(summary_state) :
1511 Well::ProductionControls(0);
1514 auto& ws = well_state_copy.well(this->index_of_well_);
1515 if (well_copy.well_ecl_.isInjector()) {
1516 inj_controls.bhp_limit = bhp;
1517 ws.injection_cmode = Well::InjectorCMode::BHP;
1519 prod_controls.bhp_limit = bhp;
1520 ws.production_cmode = Well::ProducerCMode::BHP;
1525 const int np = this->number_of_phases_;
1526 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1527 for (
int phase = 0; phase < np; ++phase){
1528 well_state_copy.wellRates(this->index_of_well_)[phase]
1529 = sign * ws.well_potentials[phase];
1531 well_copy.updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
1532 well_copy.initPrimaryVariablesEvaluation();
1533 well_copy.computeAccumWell();
1535 const double dt = simulator.timeStepSize();
1536 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1538 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1539 " potentials are computed based on unconverged solution";
1540 deferred_logger.debug(msg);
1542 well_copy.updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
1543 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1544 well_copy.initPrimaryVariablesEvaluation();
1545 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1551 template<
typename TypeTag>
1553 StandardWell<TypeTag>::
1554 computeWellPotentialWithTHP(
const Simulator& simulator,
1555 DeferredLogger& deferred_logger,
1556 const WellState &well_state)
const
1558 std::vector<double> potentials(this->number_of_phases_, 0.0);
1559 const auto& summary_state = simulator.vanguard().summaryState();
1561 const auto& well = this->well_ecl_;
1562 if (well.isInjector()){
1563 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1564 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1565 if (bhp_at_thp_limit) {
1566 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1567 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1569 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1570 "Failed in getting converged thp based potential calculation for well "
1571 + name() +
". Instead the bhp based value is used");
1572 const double bhp = controls.bhp_limit;
1573 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1576 computeWellRatesWithThpAlqProd(
1577 simulator, summary_state,
1578 deferred_logger, potentials, this->getALQ(well_state)
1585 template<
typename TypeTag>
1587 StandardWell<TypeTag>::
1588 computeWellPotentialsImplicit(
const Simulator& simulator,
1589 std::vector<double>& well_potentials,
1590 DeferredLogger& deferred_logger)
const
1595 StandardWell<TypeTag> well_copy(*
this);
1598 WellState well_state_copy = simulator.problem().wellModel().wellState();
1599 const auto& group_state = simulator.problem().wellModel().groupState();
1600 auto& ws = well_state_copy.well(this->index_of_well_);
1603 const auto& summary_state = simulator.vanguard().summaryState();
1604 auto inj_controls = well_copy.well_ecl_.isInjector()
1605 ? well_copy.well_ecl_.injectionControls(summary_state)
1606 : Well::InjectionControls(0);
1607 auto prod_controls = well_copy.well_ecl_.isProducer()
1608 ? well_copy.well_ecl_.productionControls(summary_state) :
1609 Well::ProductionControls(0);
1612 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1615 const int np = this->number_of_phases_;
1616 bool trivial =
true;
1617 for (
int phase = 0; phase < np; ++phase){
1618 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1621 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1622 for (
int phase = 0; phase < np; ++phase) {
1623 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1627 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1628 const double dt = simulator.timeStepSize();
1630 bool converged =
false;
1631 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1632 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1634 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1638 well_potentials.clear();
1639 well_potentials.resize(np, 0.0);
1640 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
1641 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
1642 well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
1648 template<
typename TypeTag>
1650 StandardWell<TypeTag>::
1651 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &simulator,
1652 const SummaryState &summary_state,
1653 DeferredLogger &deferred_logger,
1654 std::vector<double> &potentials,
1658 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1659 simulator, summary_state, alq, deferred_logger);
1660 if (bhp_at_thp_limit) {
1661 const auto& controls = this->well_ecl_.productionControls(summary_state);
1662 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1663 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1666 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1667 "Failed in getting converged thp based potential calculation for well "
1668 + name() +
". Instead the bhp based value is used");
1669 const auto& controls = this->well_ecl_.productionControls(summary_state);
1670 bhp = controls.bhp_limit;
1671 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1676 template<
typename TypeTag>
1678 StandardWell<TypeTag>::
1679 computeWellRatesWithThpAlqProd(
const Simulator& simulator,
1680 const SummaryState& summary_state,
1681 DeferredLogger& deferred_logger,
1682 std::vector<double>& potentials,
1686 computeWellRatesAndBhpWithThpAlqProd(simulator,
1693 template<
typename TypeTag>
1698 std::vector<double>& well_potentials,
1701 const auto [compute_potential, bhp_controlled_well] =
1702 this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
1704 if (!compute_potential) {
1708 bool converged_implicit =
false;
1709 if (this->param_.local_well_solver_control_switching_) {
1710 converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
1712 if (!converged_implicit) {
1714 const auto& summaryState = simulator.vanguard().summaryState();
1715 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1725 const auto& ws = well_state.well(this->index_of_well_);
1726 if (this->isInjector())
1727 bhp = std::max(ws.bhp, bhp);
1729 bhp = std::min(ws.bhp, bhp);
1731 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1732 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1735 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1739 this->checkNegativeWellPotentials(well_potentials,
1740 this->param_.check_well_operability_,
1750 template<
typename TypeTag>
1754 const int openConnIdx)
const
1756 return (openConnIdx < 0)
1758 : this->connections_.rho(openConnIdx);
1765 template<
typename TypeTag>
1767 StandardWell<TypeTag>::
1768 updatePrimaryVariables(
const SummaryState& summary_state,
1769 const WellState& well_state,
1770 DeferredLogger& deferred_logger)
1772 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1774 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
1775 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1778 if constexpr (Base::has_polymermw) {
1779 this->primary_variables_.updatePolyMW(well_state);
1782 this->primary_variables_.checkFinite(deferred_logger);
1788 template<
typename TypeTag>
1790 StandardWell<TypeTag>::
1791 getRefDensity()
const
1793 return this->connections_.rho();
1799 template<
typename TypeTag>
1801 StandardWell<TypeTag>::
1802 updateWaterMobilityWithPolymer(
const Simulator& simulator,
1804 std::vector<EvalWell>& mob,
1805 DeferredLogger& deferred_logger)
const
1807 const int cell_idx = this->well_cells_[perf];
1808 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1809 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1813 if (this->isInjector()) {
1815 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1816 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1817 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1820 if (PolymerModule::hasPlyshlog()) {
1823 if (this->isInjector() && this->wpolymer() == 0.) {
1828 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1829 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1831 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1832 PerforationRates perf_rates;
1833 double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quant, cell_idx);
1834 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1835 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1836 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1837 perf_rates, deferred_logger);
1839 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1840 const auto& material_law_manager = simulator.problem().materialLawManager();
1841 const auto& scaled_drainage_info =
1842 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1843 const double swcr = scaled_drainage_info.Swcr;
1844 const EvalWell poro = this->extendEval(int_quant.porosity());
1845 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1847 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1848 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1849 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1851 if (PolymerModule::hasShrate()) {
1854 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1856 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1857 int_quant.pvtRegionIndex(),
1860 mob[waterCompIdx] /= shear_factor;
1864 template<
typename TypeTag>
1866 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
1868 this->linSys_.extract(jacobian);
1872 template <
typename TypeTag>
1874 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1875 const BVector& weights,
1876 const int pressureVarIndex,
1877 const bool use_well_weights,
1878 const WellState& well_state)
const
1880 this->linSys_.extractCPRPressureMatrix(jacobian,
1891 template<
typename TypeTag>
1892 typename StandardWell<TypeTag>::EvalWell
1893 StandardWell<TypeTag>::
1894 pskinwater(
const double throughput,
1895 const EvalWell& water_velocity,
1896 DeferredLogger& deferred_logger)
const
1898 if constexpr (Base::has_polymermw) {
1899 const int water_table_id = this->polymerWaterTable_();
1900 if (water_table_id <= 0) {
1901 OPM_DEFLOG_THROW(std::runtime_error,
1902 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1905 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1906 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1908 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1909 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1912 OPM_DEFLOG_THROW(std::runtime_error,
1913 fmt::format(
"Polymermw is not activated, while injecting "
1914 "skin pressure is requested for well {}", name()),
1923 template<
typename TypeTag>
1924 typename StandardWell<TypeTag>::EvalWell
1925 StandardWell<TypeTag>::
1926 pskin(
const double throughput,
1927 const EvalWell& water_velocity,
1928 const EvalWell& poly_inj_conc,
1929 DeferredLogger& deferred_logger)
const
1931 if constexpr (Base::has_polymermw) {
1932 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
1933 const EvalWell water_velocity_abs = abs(water_velocity);
1934 if (poly_inj_conc == 0.) {
1935 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1937 const int polymer_table_id = this->polymerTable_();
1938 if (polymer_table_id <= 0) {
1939 OPM_DEFLOG_THROW(std::runtime_error,
1940 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1943 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1944 const double reference_concentration = skprpolytable.refConcentration;
1945 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1947 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1948 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1949 if (poly_inj_conc == reference_concentration) {
1950 return sign * pskin_poly;
1953 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1954 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1955 return sign * pskin;
1957 OPM_DEFLOG_THROW(std::runtime_error,
1958 fmt::format(
"Polymermw is not activated, while injecting "
1959 "skin pressure is requested for well {}", name()),
1968 template<
typename TypeTag>
1969 typename StandardWell<TypeTag>::EvalWell
1970 StandardWell<TypeTag>::
1971 wpolymermw(
const double throughput,
1972 const EvalWell& water_velocity,
1973 DeferredLogger& deferred_logger)
const
1975 if constexpr (Base::has_polymermw) {
1976 const int table_id = this->polymerInjTable_();
1977 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
1978 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1979 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
1980 if (this->wpolymer() == 0.) {
1981 return molecular_weight;
1983 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
1984 return molecular_weight;
1986 OPM_DEFLOG_THROW(std::runtime_error,
1987 fmt::format(
"Polymermw is not activated, while injecting "
1988 "polymer molecular weight is requested for well {}", name()),
1997 template<
typename TypeTag>
1999 StandardWell<TypeTag>::
2000 updateWaterThroughput(
const double dt, WellState &well_state)
const
2002 if constexpr (Base::has_polymermw) {
2003 if (this->isInjector()) {
2004 auto& ws = well_state.well(this->index_of_well_);
2005 auto& perf_water_throughput = ws.perf_data.water_throughput;
2006 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2007 const double perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
2009 if (perf_water_vel > 0.) {
2010 perf_water_throughput[perf] += perf_water_vel * dt;
2021 template<
typename TypeTag>
2023 StandardWell<TypeTag>::
2024 handleInjectivityRate(
const Simulator& simulator,
2026 std::vector<EvalWell>& cq_s)
const
2028 const int cell_idx = this->well_cells_[perf];
2029 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2030 const auto& fs = int_quants.fluidState();
2031 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2032 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2033 const int wat_vel_index = Bhp + 1 + perf;
2034 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2038 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2044 template<
typename TypeTag>
2046 StandardWell<TypeTag>::
2047 handleInjectivityEquations(
const Simulator& simulator,
2048 const WellState& well_state,
2050 const EvalWell& water_flux_s,
2051 DeferredLogger& deferred_logger)
2053 const int cell_idx = this->well_cells_[perf];
2054 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2055 const auto& fs = int_quants.fluidState();
2056 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2057 const EvalWell water_flux_r = water_flux_s / b_w;
2058 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2059 const EvalWell water_velocity = water_flux_r / area;
2060 const int wat_vel_index = Bhp + 1 + perf;
2063 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2065 const auto& ws = well_state.well(this->index_of_well_);
2066 const auto& perf_data = ws.perf_data;
2067 const auto& perf_water_throughput = perf_data.water_throughput;
2068 const double throughput = perf_water_throughput[perf];
2069 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2071 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2072 poly_conc.setValue(this->wpolymer());
2075 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2076 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2078 StandardWellAssemble<FluidSystem,Indices>(*this).
2079 assembleInjectivityEq(eq_pskin,
2084 this->primary_variables_.numWellEq(),
2092 template<
typename TypeTag>
2094 StandardWell<TypeTag>::
2095 checkConvergenceExtraEqs(
const std::vector<double>& res,
2096 ConvergenceReport& report)
const
2101 if constexpr (Base::has_polymermw) {
2102 WellConvergence(*this).
2103 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2111 template<
typename TypeTag>
2113 StandardWell<TypeTag>::
2114 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2115 const IntensiveQuantities& int_quants,
2116 const WellState& well_state,
2118 std::vector<RateVector>& connectionRates,
2119 DeferredLogger& deferred_logger)
const
2122 EvalWell cq_s_polymw = cq_s_poly;
2123 if (this->isInjector()) {
2124 const int wat_vel_index = Bhp + 1 + perf;
2125 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2126 if (water_velocity > 0.) {
2127 const auto& ws = well_state.well(this->index_of_well_);
2128 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2129 const double throughput = perf_water_throughput[perf];
2130 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2131 cq_s_polymw *= molecular_weight;
2137 }
else if (this->isProducer()) {
2138 if (cq_s_polymw < 0.) {
2139 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2146 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2154 template<
typename TypeTag>
2155 std::optional<double>
2156 StandardWell<TypeTag>::
2157 computeBhpAtThpLimitProd(
const WellState& well_state,
2158 const Simulator& simulator,
2159 const SummaryState& summary_state,
2160 DeferredLogger& deferred_logger)
const
2162 return computeBhpAtThpLimitProdWithAlq(simulator,
2164 this->getALQ(well_state),
2168 template<
typename TypeTag>
2169 std::optional<double>
2170 StandardWell<TypeTag>::
2171 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2172 const SummaryState& summary_state,
2173 const double alq_value,
2174 DeferredLogger& deferred_logger)
const
2177 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2183 std::vector<double> rates(3);
2184 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2185 this->adaptRatesForVFP(rates);
2189 double max_pressure = 0.0;
2190 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2191 const int cell_idx = this->well_cells_[perf];
2192 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2193 const auto& fs = int_quants.fluidState();
2194 double pressure_cell = this->getPerfCellPressure(fs).value();
2195 max_pressure = std::max(max_pressure, pressure_cell);
2197 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2200 this->connections_.rho(),
2202 this->getTHPConstraint(summary_state),
2206 auto v = frates(*bhpAtLimit);
2207 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2213 auto fratesIter = [
this, &simulator, &deferred_logger](
const double bhp) {
2217 std::vector<double> rates(3);
2218 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2219 this->adaptRatesForVFP(rates);
2223 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2226 this->connections_.rho(),
2228 this->getTHPConstraint(summary_state),
2234 auto v = frates(*bhpAtLimit);
2235 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2241 return std::nullopt;
2246 template<
typename TypeTag>
2247 std::optional<double>
2248 StandardWell<TypeTag>::
2249 computeBhpAtThpLimitInj(
const Simulator& simulator,
2250 const SummaryState& summary_state,
2251 DeferredLogger& deferred_logger)
const
2254 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2260 std::vector<double> rates(3);
2261 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2265 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2267 this->connections_.rho(),
2278 template<
typename TypeTag>
2280 StandardWell<TypeTag>::
2281 iterateWellEqWithControl(
const Simulator& simulator,
2283 const Well::InjectionControls& inj_controls,
2284 const Well::ProductionControls& prod_controls,
2285 WellState& well_state,
2286 const GroupState& group_state,
2287 DeferredLogger& deferred_logger)
2289 const int max_iter = this->param_.max_inner_iter_wells_;
2292 bool relax_convergence =
false;
2293 this->regularize_ =
false;
2294 const auto& summary_state = simulator.vanguard().summaryState();
2296 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2298 if (it > this->param_.strict_inner_iter_wells_) {
2299 relax_convergence =
true;
2300 this->regularize_ =
true;
2303 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2305 converged = report.converged();
2311 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2318 initPrimaryVariablesEvaluation();
2319 }
while (it < max_iter);
2325 template<
typename TypeTag>
2327 StandardWell<TypeTag>::
2328 iterateWellEqWithSwitching(
const Simulator& simulator,
2330 const Well::InjectionControls& inj_controls,
2331 const Well::ProductionControls& prod_controls,
2332 WellState& well_state,
2333 const GroupState& group_state,
2334 DeferredLogger& deferred_logger,
2335 const bool fixed_control ,
2336 const bool fixed_status )
2338 const int max_iter = this->param_.max_inner_iter_wells_;
2340 bool converged =
false;
2341 bool relax_convergence =
false;
2342 this->regularize_ =
false;
2343 const auto& summary_state = simulator.vanguard().summaryState();
2348 constexpr int min_its_after_switch = 4;
2349 int its_since_last_switch = min_its_after_switch;
2350 int switch_count= 0;
2352 const auto well_status_orig = this->wellStatus_;
2353 const auto operability_orig = this->operability_status_;
2354 auto well_status_cur = well_status_orig;
2355 int status_switch_count = 0;
2357 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2358 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2360 const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) &&
2361 (!fixed_control || !fixed_status) && allow_open;
2362 bool changed =
false;
2363 bool final_check =
false;
2365 this->operability_status_.resetOperability();
2366 this->operability_status_.solvable =
true;
2368 its_since_last_switch++;
2369 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2370 const double wqTotal = this->primary_variables_.eval(WQTotal).value();
2371 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger, fixed_control, fixed_status);
2373 its_since_last_switch = 0;
2375 if (well_status_cur != this->wellStatus_) {
2376 well_status_cur = this->wellStatus_;
2377 status_switch_count++;
2380 if (!changed && final_check) {
2383 final_check =
false;
2387 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2389 if (it > this->param_.strict_inner_iter_wells_) {
2390 relax_convergence =
true;
2391 this->regularize_ =
true;
2394 auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2396 converged = report.converged();
2400 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2402 its_since_last_switch = min_its_after_switch;
2409 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2410 initPrimaryVariablesEvaluation();
2412 }
while (it < max_iter);
2415 if (allow_switching){
2417 const bool is_stopped = this->wellIsStopped();
2418 if (this->wellHasTHPConstraints(summary_state)){
2419 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2420 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2422 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2426 this->wellStatus_ = well_status_orig;
2427 this->operability_status_ = operability_orig;
2428 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2429 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2430 deferred_logger.debug(message);
2436 template<
typename TypeTag>
2443 std::vector<double> well_q_s(this->num_components_, 0.);
2444 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2445 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2446 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2447 const int cell_idx = this->well_cells_[perf];
2448 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2449 std::vector<Scalar> mob(this->num_components_, 0.);
2450 getMobility(simulator, perf, mob, deferred_logger);
2451 std::vector<Scalar> cq_s(this->num_components_, 0.);
2452 double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
2453 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2454 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2456 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2457 cq_s, perf_rates, deferred_logger);
2458 for (
int comp = 0; comp < this->num_components_; ++comp) {
2459 well_q_s[comp] += cq_s[comp];
2462 const auto& comm = this->parallel_well_info_.communication();
2463 if (comm.size() > 1)
2465 comm.sum(well_q_s.data(), well_q_s.size());
2472 template <
typename TypeTag>
2477 const int num_pri_vars = this->primary_variables_.numWellEq();
2478 std::vector<double> retval(num_pri_vars);
2479 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2480 retval[ii] = this->primary_variables_.value(ii);
2489 template <
typename TypeTag>
2491 StandardWell<TypeTag>::
2492 setPrimaryVars(std::vector<double>::const_iterator it)
2494 const int num_pri_vars = this->primary_variables_.numWellEq();
2495 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2496 this->primary_variables_.setValue(ii, it[ii]);
2498 return num_pri_vars;
2502 template <
typename TypeTag>
2503 typename StandardWell<TypeTag>::Eval
2504 StandardWell<TypeTag>::
2505 connectionRateEnergy(
const double maxOilSaturation,
2506 const std::vector<EvalWell>& cq_s,
2507 const IntensiveQuantities& intQuants,
2508 DeferredLogger& deferred_logger)
const
2510 auto fs = intQuants.fluidState();
2512 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2513 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2518 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2519 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2520 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2521 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2522 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2525 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2526 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2531 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2533 deferred_logger.debug(
2534 fmt::format(
"Problematic d value {} obtained for well {}"
2535 " during calculateSinglePerf with rs {}"
2536 ", rv {}. Continue as if no dissolution (rs = 0) and"
2537 " vaporization (rv = 0) for this connection.",
2538 d, this->name(), fs.Rs(), fs.Rv()));
2539 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2541 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2542 cq_r_thermal = (cq_s[gasCompIdx] -
2543 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2544 (d * this->extendEval(fs.invB(phaseIdx)) );
2545 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2547 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2549 (d * this->extendEval(fs.invB(phaseIdx)) );
2555 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
2557 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2558 fs.setTemperature(this->well_ecl_.temperature());
2559 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2560 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2561 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2562 paramCache.setRegionIndex(pvtRegionIdx);
2563 paramCache.setMaxOilSat(maxOilSaturation);
2564 paramCache.updatePhase(fs, phaseIdx);
2566 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2567 fs.setDensity(phaseIdx, rho);
2568 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2569 fs.setEnthalpy(phaseIdx, h);
2570 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2571 result += getValue(cq_r_thermal);
2574 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2575 result += Base::restrictEval(cq_r_thermal);
2583 template <
typename TypeTag>
2584 template<
class Value>
2586 StandardWell<TypeTag>::
2587 gasOilPerfRateInj(
const std::vector<Value>& cq_s,
2588 PerforationRates& perf_rates,
2591 const Value& pressure,
2593 DeferredLogger& deferred_logger)
const
2595 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2596 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2605 const double d = 1.0 - getValue(rv) * getValue(rs);
2608 deferred_logger.debug(dValueError(d, this->name(),
2609 "gasOilPerfRateInj",
2614 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2617 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
2620 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2625 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2631 template <
typename TypeTag>
2632 template<
class Value>
2634 StandardWell<TypeTag>::
2635 gasOilPerfRateProd(std::vector<Value>& cq_s,
2636 PerforationRates& perf_rates,
2639 const Value& rvw)
const
2641 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2642 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2643 const Value cq_sOil = cq_s[oilCompIdx];
2644 const Value cq_sGas = cq_s[gasCompIdx];
2645 const Value dis_gas = rs * cq_sOil;
2646 const Value vap_oil = rv * cq_sGas;
2648 cq_s[gasCompIdx] += dis_gas;
2649 cq_s[oilCompIdx] += vap_oil;
2652 if (this->isProducer()) {
2653 perf_rates.dis_gas = getValue(dis_gas);
2654 perf_rates.vap_oil = getValue(vap_oil);
2657 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2658 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2659 const Value vap_wat = rvw * cq_sGas;
2660 cq_s[waterCompIdx] += vap_wat;
2661 if (this->isProducer())
2662 perf_rates.vap_wat = getValue(vap_wat);
2667 template <
typename TypeTag>
2668 template<
class Value>
2670 StandardWell<TypeTag>::
2671 gasWaterPerfRateProd(std::vector<Value>& cq_s,
2672 PerforationRates& perf_rates,
2674 const Value& rsw)
const
2676 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2677 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2678 const Value cq_sWat = cq_s[waterCompIdx];
2679 const Value cq_sGas = cq_s[gasCompIdx];
2680 const Value vap_wat = rvw * cq_sGas;
2681 const Value dis_gas_wat = rsw * cq_sWat;
2682 cq_s[waterCompIdx] += vap_wat;
2683 cq_s[gasCompIdx] += dis_gas_wat;
2684 if (this->isProducer()) {
2685 perf_rates.vap_wat = getValue(vap_wat);
2686 perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
2691 template <
typename TypeTag>
2692 template<
class Value>
2694 StandardWell<TypeTag>::
2695 gasWaterPerfRateInj(
const std::vector<Value>& cq_s,
2696 PerforationRates& perf_rates,
2699 const Value& pressure,
2700 DeferredLogger& deferred_logger)
const
2703 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2704 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2706 const double dw = 1.0 - getValue(rvw) * getValue(rsw);
2709 deferred_logger.debug(dValueError(dw, this->name(),
2710 "gasWaterPerfRateInj",
2711 rsw, rvw, pressure));
2715 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
2718 perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
2723 template <
typename TypeTag>
2724 template<
class Value>
2726 StandardWell<TypeTag>::
2727 disOilVapWatVolumeRatio(Value& volumeRatio,
2730 const Value& pressure,
2731 const std::vector<Value>& cmix_s,
2732 const std::vector<Value>& b_perfcells_dense,
2733 DeferredLogger& deferred_logger)
const
2735 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2736 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2738 const Value d = 1.0 - rvw * rsw;
2741 deferred_logger.debug(dValueError(d, this->name(),
2742 "disOilVapWatVolumeRatio",
2743 rsw, rvw, pressure));
2745 const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
2746 : cmix_s[waterCompIdx];
2747 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
2749 const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
2750 : cmix_s[waterCompIdx];
2751 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
2755 template <
typename TypeTag>
2756 template<
class Value>
2758 StandardWell<TypeTag>::
2759 gasOilVolumeRatio(Value& volumeRatio,
2762 const Value& pressure,
2763 const std::vector<Value>& cmix_s,
2764 const std::vector<Value>& b_perfcells_dense,
2765 DeferredLogger& deferred_logger)
const
2767 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2768 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2770 const Value d = 1.0 - rv * rs;
2773 deferred_logger.debug(dValueError(d, this->name(),
2774 "gasOilVolumeRatio",
2777 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
2778 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
2780 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
2781 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
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:59
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition StandardWell_impl.hpp:1391
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:61
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Definition PerforationData.hpp:40