22#ifndef OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
23#define OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
26#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
28#include <opm/simulators/wells/StandardWell.hpp>
31#include <opm/common/Exceptions.hpp>
33#include <opm/input/eclipse/Units/Units.hpp>
35#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
36#include <opm/simulators/wells/StandardWellAssemble.hpp>
37#include <opm/simulators/wells/VFPHelpers.hpp>
38#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
39#include <opm/simulators/wells/WellConvergence.hpp>
45#include <fmt/format.h>
50 template<
typename TypeTag>
51 StandardWell<TypeTag>::
52 StandardWell(
const Well& well,
55 const ModelParameters& param,
57 const int pvtRegionIdx,
66 assert(this->num_components_ == numWellConservationEq);
73 template<
typename TypeTag>
75 StandardWell<TypeTag>::
76 init(
const PhaseUsage* phase_usage_arg,
77 const std::vector<Scalar>& depth_arg,
78 const Scalar gravity_arg,
79 const std::vector< Scalar >& B_avg,
80 const bool changed_to_open_this_step)
82 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
83 this->StdWellEval::init(this->perf_depth_, depth_arg, Base::has_polymermw);
90 template<
typename TypeTag>
93 StandardWell<TypeTag>::
94 computePerfRate(
const IntensiveQuantities& intQuants,
95 const std::vector<Value>& mob,
97 const std::vector<Scalar>& Tw,
100 std::vector<Value>& cq_s,
101 PerforationRates<Scalar>& perf_rates,
102 DeferredLogger& deferred_logger)
const
104 auto obtain = [
this](
const Eval& value)
106 if constexpr (std::is_same_v<Value, Scalar>) {
107 static_cast<void>(
this);
108 return getValue(value);
110 return this->extendEval(value);
113 auto obtainN = [](
const auto& value)
115 if constexpr (std::is_same_v<Value, Scalar>) {
116 return getValue(value);
121 auto zeroElem = [
this]()
123 if constexpr (std::is_same_v<Value, Scalar>) {
124 static_cast<void>(
this);
127 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
131 const auto& fs = intQuants.fluidState();
132 const Value pressure = obtain(this->getPerfCellPressure(fs));
133 const Value rs = obtain(fs.Rs());
134 const Value rv = obtain(fs.Rv());
135 const Value rvw = obtain(fs.Rvw());
136 const Value rsw = obtain(fs.Rsw());
138 std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
139 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
140 if (!FluidSystem::phaseIsActive(phaseIdx)) {
143 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
144 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
146 if constexpr (has_solvent) {
147 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
150 if constexpr (has_zFraction) {
151 if (this->isInjector()) {
152 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
153 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
154 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
158 Value skin_pressure = zeroElem();
160 if (this->isInjector()) {
161 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
162 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
167 std::vector<Value> cmix_s(this->numComponents(), zeroElem());
168 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
169 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
190 template<
typename TypeTag>
191 template<
class Value>
193 StandardWell<TypeTag>::
194 computePerfRate(
const std::vector<Value>& mob,
195 const Value& pressure,
201 std::vector<Value>& b_perfcells_dense,
202 const std::vector<Scalar>& Tw,
205 const Value& skin_pressure,
206 const std::vector<Value>& cmix_s,
207 std::vector<Value>& cq_s,
208 PerforationRates<Scalar>& perf_rates,
209 DeferredLogger& deferred_logger)
const
212 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
213 Value drawdown = pressure - well_pressure;
214 if (this->isInjector()) {
215 drawdown += skin_pressure;
218 RatioCalculator<Value> ratioCalc{
219 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
220 ? Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)
222 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
223 ? Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)
225 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
226 ? Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)
234 if (!allow_cf && this->isInjector()) {
239 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
240 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
241 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
244 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
245 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
247 ratioCalc.gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw,
248 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
250 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
251 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
253 ratioCalc.gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw, this->isProducer());
257 if (!allow_cf && this->isProducer()) {
262 Value total_mob_dense = mob[0];
263 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
264 total_mob_dense += mob[componentIdx];
268 Value volumeRatio = bhp * 0.0;
270 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
271 ratioCalc.disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
272 cmix_s, b_perfcells_dense, deferred_logger);
276 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
277 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
278 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
281 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
282 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
283 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) &&
291 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
293 ratioCalc.gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
294 cmix_s, b_perfcells_dense,
297 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
298 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
299 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
301 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
302 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
303 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
309 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
310 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
311 Value cqt_is = cqt_i / volumeRatio;
312 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
316 if (this->isProducer()) {
317 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
318 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
320 ratioCalc.gasOilPerfRateInj(cq_s, perf_rates,
321 rv, rs, pressure, rvw,
322 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
325 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
326 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
329 ratioCalc.gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
330 pressure, deferred_logger);
337 template<
typename TypeTag>
339 StandardWell<TypeTag>::
340 assembleWellEqWithoutIteration(
const Simulator& simulator,
342 const Well::InjectionControls& inj_controls,
343 const Well::ProductionControls& prod_controls,
344 WellState<Scalar>& well_state,
345 const GroupState<Scalar>& group_state,
346 DeferredLogger& deferred_logger)
350 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
353 this->linSys_.clear();
355 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
356 prod_controls, well_state,
357 group_state, deferred_logger);
363 template<
typename TypeTag>
365 StandardWell<TypeTag>::
366 assembleWellEqWithoutIterationImpl(
const Simulator& simulator,
368 const Well::InjectionControls& inj_controls,
369 const Well::ProductionControls& prod_controls,
370 WellState<Scalar>& well_state,
371 const GroupState<Scalar>& group_state,
372 DeferredLogger& deferred_logger)
375 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
376 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
378 auto& ws = well_state.well(this->index_of_well_);
379 ws.phase_mixing_rates.fill(0.0);
382 const int np = this->number_of_phases_;
384 std::vector<RateVector> connectionRates = this->connectionRates_;
386 auto& perf_data = ws.perf_data;
387 auto& perf_rates = perf_data.phase_rates;
388 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
390 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
391 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
392 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
393 calculateSinglePerf(simulator, perf, well_state, connectionRates,
394 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
397 if constexpr (has_polymer && Base::has_polymermw) {
398 if (this->isInjector()) {
399 handleInjectivityEquations(simulator, well_state, perf,
400 water_flux_s, deferred_logger);
403 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
405 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
407 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
409 StandardWellAssemble<FluidSystem,Indices>(*this).
410 assemblePerforationEq(cq_s_effective,
413 this->primary_variables_.numWellEq(),
417 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
418 auto& perf_rate_solvent = perf_data.solvent_rates;
419 perf_rate_solvent[perf] = cq_s[componentIdx].value();
421 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
425 if constexpr (has_zFraction) {
426 StandardWellAssemble<FluidSystem,Indices>(*this).
427 assembleZFracEq(cq_s_zfrac_effective,
429 this->primary_variables_.numWellEq(),
434 this->connectionRates_ = connectionRates;
439 const auto& comm = this->parallel_well_info_.communication();
440 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
444 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
447 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
450 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
451 if (FluidSystem::numActivePhases() > 1) {
453 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
454 this->F0_[componentIdx]) * volume / dt;
456 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
457 StandardWellAssemble<FluidSystem,Indices>(*this).
458 assembleSourceEq(resWell_loc,
460 this->primary_variables_.numWellEq(),
464 const auto& summaryState = simulator.vanguard().summaryState();
465 const Schedule& schedule = simulator.vanguard().schedule();
466 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
467 StandardWellAssemble<FluidSystem,Indices>(*this).
468 assembleControlEq(well_state, group_state,
469 schedule, summaryState,
470 inj_controls, prod_controls,
471 this->primary_variables_,
472 this->connections_.rho(),
474 stopped_or_zero_target,
480 this->linSys_.invert();
482 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
489 template<
typename TypeTag>
491 StandardWell<TypeTag>::
492 calculateSinglePerf(
const Simulator& simulator,
494 WellState<Scalar>& well_state,
495 std::vector<RateVector>& connectionRates,
496 std::vector<EvalWell>& cq_s,
497 EvalWell& water_flux_s,
498 EvalWell& cq_s_zfrac_effective,
499 DeferredLogger& deferred_logger)
const
501 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
502 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
503 const int cell_idx = this->well_cells_[perf];
504 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
505 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
506 getMobility(simulator, perf, mob, deferred_logger);
508 PerforationRates<Scalar> perf_rates;
509 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
510 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
511 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
512 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
513 cq_s, perf_rates, deferred_logger);
515 auto& ws = well_state.well(this->index_of_well_);
516 auto& perf_data = ws.perf_data;
517 if constexpr (has_polymer && Base::has_polymermw) {
518 if (this->isInjector()) {
521 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
522 water_flux_s = cq_s[water_comp_idx];
525 handleInjectivityRate(simulator, perf, cq_s);
530 if (this->isProducer()) {
531 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
532 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
533 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
534 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
535 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.dis_gas;
536 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.dis_gas_in_water;
537 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.vap_oil;
538 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.vap_wat;
541 if constexpr (has_energy) {
542 connectionRates[perf][Indices::contiEnergyEqIdx] =
543 connectionRateEnergy(simulator.problem().maxOilSaturation(cell_idx),
544 cq_s, intQuants, deferred_logger);
547 if constexpr (has_polymer) {
548 std::variant<Scalar,EvalWell> polymerConcentration;
549 if (this->isInjector()) {
550 polymerConcentration = this->wpolymer();
552 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
553 intQuants.polymerViscosityCorrection());
556 [[maybe_unused]] EvalWell cq_s_poly;
557 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
559 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
560 cq_s, polymerConcentration);
562 if constexpr (Base::has_polymermw) {
563 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
564 perf, connectionRates, deferred_logger);
568 if constexpr (has_foam) {
569 std::variant<Scalar,EvalWell> foamConcentration;
570 if (this->isInjector()) {
571 foamConcentration = this->wfoam();
573 foamConcentration = this->extendEval(intQuants.foamConcentration());
575 connectionRates[perf][Indices::contiFoamEqIdx] =
576 this->connections_.connectionRateFoam(cq_s, foamConcentration,
577 FoamModule::transportPhase(),
581 if constexpr (has_zFraction) {
582 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
583 if (this->isInjector()) {
584 solventConcentration = this->wsolvent();
586 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
587 this->extendEval(intQuants.yVolume())};
589 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
590 cq_s_zfrac_effective) =
591 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
592 perf_rates.dis_gas, cq_s,
593 solventConcentration);
596 if constexpr (has_brine) {
597 std::variant<Scalar,EvalWell> saltConcentration;
598 if (this->isInjector()) {
599 saltConcentration = this->wsalt();
601 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
604 connectionRates[perf][Indices::contiBrineEqIdx] =
605 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
606 perf_rates.vap_wat, cq_s,
610 if constexpr (has_micp) {
611 std::variant<Scalar,EvalWell> microbialConcentration;
612 std::variant<Scalar,EvalWell> oxygenConcentration;
613 std::variant<Scalar,EvalWell> ureaConcentration;
614 if (this->isInjector()) {
615 microbialConcentration = this->wmicrobes();
616 oxygenConcentration = this->woxygen();
617 ureaConcentration = this->wurea();
619 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
620 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
621 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
623 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
624 connectionRates[perf][Indices::contiOxygenEqIdx],
625 connectionRates[perf][Indices::contiUreaEqIdx]) =
626 this->connections_.connectionRatesMICP(cq_s,
627 microbialConcentration,
633 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
638 template<
typename TypeTag>
639 template<
class Value>
641 StandardWell<TypeTag>::
642 getMobility(
const Simulator& simulator,
644 std::vector<Value>& mob,
645 DeferredLogger& deferred_logger)
const
647 auto obtain = [
this](
const Eval& value)
649 if constexpr (std::is_same_v<Value, Scalar>) {
650 static_cast<void>(
this);
651 return getValue(value);
653 return this->extendEval(value);
656 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
657 obtain, deferred_logger);
660 if constexpr (has_polymer) {
661 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
662 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
667 if constexpr (!Base::has_polymermw) {
668 if constexpr (std::is_same_v<Value, Scalar>) {
669 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
670 for (std::size_t i = 0; i < mob.size(); ++i) {
671 mob_eval[i].setValue(mob[i]);
673 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
674 for (std::size_t i = 0; i < mob.size(); ++i) {
675 mob[i] = getValue(mob_eval[i]);
678 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
684 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
685 const Scalar bhp = this->primary_variables_.value(Bhp);
686 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
687 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
688 for (std::size_t i = 0; i < mob.size(); ++i) {
689 mob[i] *= multiplier;
695 template<
typename TypeTag>
697 StandardWell<TypeTag>::
698 updateWellState(
const Simulator& simulator,
699 const BVectorWell& dwells,
700 WellState<Scalar>& well_state,
701 DeferredLogger& deferred_logger)
703 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
705 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
706 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
708 const auto& summary_state = simulator.vanguard().summaryState();
709 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
710 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
717 template<
typename TypeTag>
719 StandardWell<TypeTag>::
720 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
721 const bool stop_or_zero_rate_target,
722 DeferredLogger& deferred_logger)
724 const Scalar dFLimit = this->param_.dwell_fraction_max_;
725 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
726 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
729 if constexpr (Base::has_polymermw) {
730 this->primary_variables_.updateNewtonPolyMW(dwells);
733 this->primary_variables_.checkFinite(deferred_logger);
740 template<
typename TypeTag>
742 StandardWell<TypeTag>::
743 updateWellStateFromPrimaryVariables(WellState<Scalar>& well_state,
744 const SummaryState& summary_state,
745 DeferredLogger& deferred_logger)
const
747 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
750 if constexpr (Base::has_polymermw) {
751 this->primary_variables_.copyToWellStatePolyMW(well_state);
759 template<
typename TypeTag>
761 StandardWell<TypeTag>::
762 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
767 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
768 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
770 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
771 std::vector<Scalar> mob(this->num_components_, 0.0);
772 getMobility(simulator, perf, mob, deferred_logger);
774 const int cell_idx = this->well_cells_[perf];
775 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
776 const auto& fs = int_quantities.fluidState();
778 Scalar p_r = this->getPerfCellPressure(fs).value();
781 std::vector<Scalar> b_perf(this->num_components_);
782 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
783 if (!FluidSystem::phaseIsActive(phase)) {
786 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
787 b_perf[comp_idx] = fs.invB(phase).value();
789 if constexpr (has_solvent) {
790 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
794 const Scalar h_perf = this->connections_.pressure_diff(perf);
795 const Scalar pressure_diff = p_r - h_perf;
800 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
801 deferred_logger.debug(
"CROSSFLOW_IPR",
802 "cross flow found when updateIPR for well " + name()
803 +
" . The connection is ignored in IPR calculations");
809 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
810 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
811 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
815 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
816 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
817 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
818 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
819 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
820 ipr_b_perf[comp_idx] += tw_mob;
824 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
825 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
826 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
827 const Scalar rs = (fs.Rs()).value();
828 const Scalar rv = (fs.Rv()).value();
830 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
831 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
833 ipr_a_perf[gas_comp_idx] += dis_gas_a;
834 ipr_a_perf[oil_comp_idx] += vap_oil_a;
836 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
837 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
839 ipr_b_perf[gas_comp_idx] += dis_gas_b;
840 ipr_b_perf[oil_comp_idx] += vap_oil_b;
843 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
844 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
845 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
848 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
849 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
852 template<
typename TypeTag>
854 StandardWell<TypeTag>::
855 updateIPRImplicit(
const Simulator& simulator,
856 WellState<Scalar>& well_state,
857 DeferredLogger& deferred_logger)
866 auto rates = well_state.well(this->index_of_well_).surface_rates;
868 for (std::size_t p = 0; p < rates.size(); ++p) {
869 zero_rates &= rates[p] == 0.0;
871 auto& ws = well_state.well(this->index_of_well_);
873 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
874 deferred_logger.debug(msg);
886 const auto& group_state = simulator.problem().wellModel().groupState();
888 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
889 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
891 auto inj_controls = Well::InjectionControls(0);
892 auto prod_controls = Well::ProductionControls(0);
893 prod_controls.addControl(Well::ProducerCMode::BHP);
894 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
897 const auto cmode = ws.production_cmode;
898 ws.production_cmode = Well::ProducerCMode::BHP;
899 const double dt = simulator.timeStepSize();
900 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
902 const size_t nEq = this->primary_variables_.numWellEq();
906 for (
size_t i=0; i < nEq; ++i){
911 BVectorWell x_well(1);
912 x_well[0].resize(nEq);
913 this->linSys_.solve(rhs, x_well);
915 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
916 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
917 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
918 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
920 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
922 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
925 ws.production_cmode = cmode;
928 template<
typename TypeTag>
930 StandardWell<TypeTag>::
931 checkOperabilityUnderBHPLimit(
const WellState<Scalar>& well_state,
932 const Simulator& simulator,
933 DeferredLogger& deferred_logger)
935 const auto& summaryState = simulator.vanguard().summaryState();
936 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
939 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
940 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
943 Scalar total_ipr_mass_rate = 0.0;
944 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
946 if (!FluidSystem::phaseIsActive(phaseIdx)) {
950 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
951 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
953 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
954 total_ipr_mass_rate += ipr_rate * rho;
956 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
957 this->operability_status_.operable_under_only_bhp_limit =
false;
961 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
965 std::vector<Scalar> well_rates_bhp_limit;
966 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
968 this->adaptRatesForVFP(well_rates_bhp_limit);
969 const Scalar thp_limit = this->getTHPConstraint(summaryState);
970 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
972 this->connections_.rho(),
973 this->getALQ(well_state),
976 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
977 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
988 this->operability_status_.operable_under_only_bhp_limit =
true;
989 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
997 template<
typename TypeTag>
999 StandardWell<TypeTag>::
1000 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
1001 const WellState<Scalar>& well_state,
1002 DeferredLogger& deferred_logger)
1004 const auto& summaryState = simulator.vanguard().summaryState();
1005 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1006 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1009 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1011 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1012 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1013 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1015 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1016 if (this->isProducer() && *obtain_bhp < thp_limit) {
1017 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1018 +
" bars is SMALLER than thp limit "
1019 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1020 +
" bars as a producer for well " + name();
1021 deferred_logger.debug(msg);
1023 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1024 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1025 +
" bars is LARGER than thp limit "
1026 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1027 +
" bars as a injector for well " + name();
1028 deferred_logger.debug(msg);
1031 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1032 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1033 if (!this->wellIsStopped()) {
1034 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1035 deferred_logger.debug(
" could not find bhp value at thp limit "
1036 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1037 +
" bar for well " + name() +
", the well might need to be closed ");
1046 template<
typename TypeTag>
1048 StandardWell<TypeTag>::
1049 allDrawDownWrongDirection(
const Simulator& simulator)
const
1051 bool all_drawdown_wrong_direction =
true;
1053 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1054 const int cell_idx = this->well_cells_[perf];
1055 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1056 const auto& fs = intQuants.fluidState();
1058 const Scalar pressure = this->getPerfCellPressure(fs).value();
1059 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1062 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1063 const Scalar drawdown = pressure - well_pressure;
1068 if ( (drawdown < 0. && this->isInjector()) ||
1069 (drawdown > 0. && this->isProducer()) ) {
1070 all_drawdown_wrong_direction =
false;
1075 const auto& comm = this->parallel_well_info_.communication();
1076 if (comm.size() > 1)
1078 all_drawdown_wrong_direction =
1079 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1082 return all_drawdown_wrong_direction;
1088 template<
typename TypeTag>
1090 StandardWell<TypeTag>::
1091 canProduceInjectWithCurrentBhp(
const Simulator& simulator,
1092 const WellState<Scalar>& well_state,
1093 DeferredLogger& deferred_logger)
1095 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
1096 std::vector<Scalar> well_rates;
1097 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1099 const Scalar sign = (this->isProducer()) ? -1. : 1.;
1100 const Scalar threshold = sign * std::numeric_limits<Scalar>::min();
1102 bool can_produce_inject =
false;
1103 for (
const auto value : well_rates) {
1104 if (this->isProducer() && value < threshold) {
1105 can_produce_inject =
true;
1107 }
else if (this->isInjector() && value > threshold) {
1108 can_produce_inject =
true;
1113 if (!can_produce_inject) {
1114 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1117 return can_produce_inject;
1124 template<
typename TypeTag>
1126 StandardWell<TypeTag>::
1127 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1129 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1135 template<
typename TypeTag>
1136 typename StandardWell<TypeTag>::WellConnectionProps
1137 StandardWell<TypeTag>::
1138 computePropertiesForWellConnectionPressures(
const Simulator& simulator,
1139 const WellState<Scalar>& well_state)
const
1141 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1143 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1145 return model.intensiveQuantities(cell_idx, 0)
1146 .fluidState().temperature(phase_idx).value();
1150 [&model = simulator.model()](
int cell_idx)
1152 return model.intensiveQuantities(cell_idx, 0)
1153 .fluidState().saltConcentration().value();
1157 [&model = simulator.model()](
int cell_idx)
1159 return model.intensiveQuantities(cell_idx, 0)
1160 .fluidState().pvtRegionIndex();
1164 if constexpr (Indices::enableSolvent) {
1165 prop_func.solventInverseFormationVolumeFactor =
1166 [&model = simulator.model()](
int cell_idx)
1168 return model.intensiveQuantities(cell_idx, 0)
1169 .solventInverseFormationVolumeFactor().value();
1172 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1174 return model.intensiveQuantities(cell_idx, 0)
1175 .solventRefDensity();
1179 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1186 template<
typename TypeTag>
1191 const std::vector<Scalar>& B_avg,
1193 const bool relax_tolerance)
const
1197 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1199 Scalar tol_wells = this->param_.tolerance_wells_;
1201 constexpr Scalar stopped_factor = 1.e-4;
1203 constexpr Scalar dynamic_thp_factor = 1.e-1;
1204 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1205 tol_wells = tol_wells*stopped_factor;
1206 }
else if (this->getDynamicThpLimit()) {
1207 tol_wells = tol_wells*dynamic_thp_factor;
1210 std::vector<Scalar> res;
1213 this->param_.max_residual_allowed_,
1215 this->param_.relaxed_tolerance_flow_well_,
1217 this->wellIsStopped(),
1221 checkConvergenceExtraEqs(res, report);
1230 template<
typename TypeTag>
1238 auto fluidState = [&simulator,
this](
const int perf)
1240 const auto cell_idx = this->well_cells_[perf];
1241 return simulator.
model()
1242 .intensiveQuantities(cell_idx, 0).fluidState();
1245 const int np = this->number_of_phases_;
1246 auto setToZero = [np](Scalar* x) ->
void
1248 std::fill_n(x, np, 0.0);
1251 auto addVector = [np](
const Scalar* src, Scalar* dest) ->
void
1253 std::transform(src, src + np, dest, dest, std::plus<>{});
1256 auto& ws = well_state.well(this->index_of_well_);
1257 auto& perf_data = ws.perf_data;
1258 auto* wellPI = ws.productivity_index.data();
1259 auto* connPI = perf_data.prod_index.data();
1263 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1264 auto subsetPerfID = 0;
1266 for (
const auto& perf : *this->perf_data_) {
1267 auto allPerfID = perf.ecl_index;
1269 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) -> Scalar
1274 std::vector<Scalar> mob(this->num_components_, 0.0);
1275 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1277 const auto& fs = fluidState(subsetPerfID);
1280 if (this->isInjector()) {
1281 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1282 mob, connPI, deferred_logger);
1285 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1288 addVector(connPI, wellPI);
1295 const auto& comm = this->parallel_well_info_.communication();
1296 if (comm.size() > 1) {
1297 comm.sum(wellPI, np);
1300 assert ((
static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1301 "Internal logic error in processing connections for PI/II");
1306 template<
typename TypeTag>
1307 void StandardWell<TypeTag>::
1308 computeWellConnectionDensitesPressures(
const Simulator& simulator,
1309 const WellState<Scalar>& well_state,
1310 const WellConnectionProps& props,
1311 DeferredLogger& deferred_logger)
1316 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1321 [&model = simulator.model()](
const int cell,
1322 const std::vector<int>& phases,
1323 std::vector<Scalar>& mob)
1325 const auto& iq = model.intensiveQuantities(cell, 0);
1327 std::transform(phases.begin(), phases.end(), mob.begin(),
1328 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1333 [&model = simulator.model()](
const int cell,
1334 const std::vector<int>& phases,
1335 std::vector<Scalar>& rho)
1337 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1339 std::transform(phases.begin(), phases.end(), rho.begin(),
1340 [&fs](
const int phase) { return fs.density(phase).value(); });
1344 const auto stopped_or_zero_rate_target = this->
1345 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1348 .computeProperties(stopped_or_zero_rate_target, well_state,
1349 prop_func, props, deferred_logger);
1356 template<
typename TypeTag>
1358 StandardWell<TypeTag>::
1359 computeWellConnectionPressures(
const Simulator& simulator,
1360 const WellState<Scalar>& well_state,
1361 DeferredLogger& deferred_logger)
1363 const auto props = computePropertiesForWellConnectionPressures
1364 (simulator, well_state);
1366 computeWellConnectionDensitesPressures(simulator, well_state,
1367 props, deferred_logger);
1374 template<
typename TypeTag>
1376 StandardWell<TypeTag>::
1377 solveEqAndUpdateWellState(
const Simulator& simulator,
1378 WellState<Scalar>& well_state,
1379 DeferredLogger& deferred_logger)
1381 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1385 BVectorWell dx_well(1);
1386 dx_well[0].resize(this->primary_variables_.numWellEq());
1387 this->linSys_.solve( dx_well);
1389 updateWellState(simulator, dx_well, well_state, deferred_logger);
1396 template<
typename TypeTag>
1398 StandardWell<TypeTag>::
1399 calculateExplicitQuantities(
const Simulator& simulator,
1400 const WellState<Scalar>& well_state,
1401 DeferredLogger& deferred_logger)
1403 updatePrimaryVariables(simulator, well_state, deferred_logger);
1404 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1405 this->computeAccumWell();
1410 template<
typename TypeTag>
1413 apply(
const BVector& x, BVector& Ax)
const
1415 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1417 if (this->param_.matrix_add_well_contributions_)
1423 this->linSys_.
apply(x, Ax);
1429 template<
typename TypeTag>
1432 apply(BVector& r)
const
1434 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1436 this->linSys_.
apply(r);
1442 template<
typename TypeTag>
1450 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1453 xw[0].resize(this->primary_variables_.numWellEq());
1455 this->linSys_.recoverSolutionWell(x, xw);
1456 updateWellState(simulator, xw, well_state, deferred_logger);
1462 template<
typename TypeTag>
1467 std::vector<Scalar>& well_flux,
1471 const int np = this->number_of_phases_;
1472 well_flux.resize(np, 0.0);
1474 const bool allow_cf = this->getAllowCrossFlow();
1476 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1477 const int cell_idx = this->well_cells_[perf];
1478 const auto& intQuants = simulator.
model().intensiveQuantities(cell_idx, 0);
1480 std::vector<Scalar> mob(this->num_components_, 0.);
1481 getMobility(simulator, perf, mob, deferred_logger);
1482 Scalar trans_mult = simulator.
problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1483 const auto& wellstate_nupcol = simulator.
problem().wellModel().nupcolWellState().well(this->index_of_well_);
1484 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1486 std::vector<Scalar> cq_s(this->num_components_, 0.);
1488 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1489 cq_s, perf_rates, deferred_logger);
1491 for(
int p = 0; p < np; ++p) {
1492 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1496 if constexpr (has_solvent) {
1497 const auto& pu = this->phaseUsage();
1498 assert(pu.phase_used[Gas]);
1499 const int gas_pos = pu.phase_pos[Gas];
1500 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1503 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1508 template<
typename TypeTag>
1510 StandardWell<TypeTag>::
1511 computeWellRatesWithBhpIterations(
const Simulator& simulator,
1513 std::vector<Scalar>& well_flux,
1514 DeferredLogger& deferred_logger)
const
1518 StandardWell<TypeTag> well_copy(*
this);
1519 well_copy.resetDampening();
1524 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1525 const auto& group_state = simulator.problem().wellModel().groupState();
1528 const auto& summary_state = simulator.vanguard().summaryState();
1529 auto inj_controls = well_copy.well_ecl_.isInjector()
1530 ? well_copy.well_ecl_.injectionControls(summary_state)
1531 : Well::InjectionControls(0);
1532 auto prod_controls = well_copy.well_ecl_.isProducer()
1533 ? well_copy.well_ecl_.productionControls(summary_state) :
1534 Well::ProductionControls(0);
1537 auto& ws = well_state_copy.well(this->index_of_well_);
1538 if (well_copy.well_ecl_.isInjector()) {
1539 inj_controls.bhp_limit = bhp;
1540 ws.injection_cmode = Well::InjectorCMode::BHP;
1542 prod_controls.bhp_limit = bhp;
1543 ws.production_cmode = Well::ProducerCMode::BHP;
1548 const int np = this->number_of_phases_;
1549 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1550 for (
int phase = 0; phase < np; ++phase){
1551 well_state_copy.wellRates(this->index_of_well_)[phase]
1552 = sign * ws.well_potentials[phase];
1554 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1555 well_copy.computeAccumWell();
1557 const double dt = simulator.timeStepSize();
1558 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1560 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1561 " potentials are computed based on unconverged solution";
1562 deferred_logger.debug(msg);
1564 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1565 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1566 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1572 template<
typename TypeTag>
1573 std::vector<typename StandardWell<TypeTag>::Scalar>
1574 StandardWell<TypeTag>::
1575 computeWellPotentialWithTHP(
const Simulator& simulator,
1576 DeferredLogger& deferred_logger,
1577 const WellState<Scalar>& well_state)
const
1579 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1580 const auto& summary_state = simulator.vanguard().summaryState();
1582 const auto& well = this->well_ecl_;
1583 if (well.isInjector()){
1584 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1585 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1586 if (bhp_at_thp_limit) {
1587 const Scalar bhp = std::min(*bhp_at_thp_limit,
1588 static_cast<Scalar
>(controls.bhp_limit));
1589 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1591 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1592 "Failed in getting converged thp based potential calculation for well "
1593 + name() +
". Instead the bhp based value is used");
1594 const Scalar bhp = controls.bhp_limit;
1595 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1598 computeWellRatesWithThpAlqProd(
1599 simulator, summary_state,
1600 deferred_logger, potentials, this->getALQ(well_state)
1607 template<
typename TypeTag>
1609 StandardWell<TypeTag>::
1610 computeWellPotentialsImplicit(
const Simulator& simulator,
1611 const WellState<Scalar>& well_state,
1612 std::vector<Scalar>& well_potentials,
1613 DeferredLogger& deferred_logger)
const
1618 StandardWell<TypeTag> well_copy(*
this);
1621 WellState<Scalar> well_state_copy = well_state;
1622 const auto& group_state = simulator.problem().wellModel().groupState();
1623 auto& ws = well_state_copy.well(this->index_of_well_);
1626 const auto& summary_state = simulator.vanguard().summaryState();
1627 auto inj_controls = well_copy.well_ecl_.isInjector()
1628 ? well_copy.well_ecl_.injectionControls(summary_state)
1629 : Well::InjectionControls(0);
1630 auto prod_controls = well_copy.well_ecl_.isProducer()
1631 ? well_copy.well_ecl_.productionControls(summary_state) :
1632 Well::ProductionControls(0);
1635 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1638 const int num_perf = ws.perf_data.size();
1639 for (
int perf = 0; perf < num_perf; ++perf) {
1640 ws.perf_data.pressure[perf] = ws.bhp + well_copy.connections_.pressure_diff(perf);
1643 const int np = this->number_of_phases_;
1644 bool trivial =
true;
1645 for (
int phase = 0; phase < np; ++phase){
1646 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1649 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1650 for (
int phase = 0; phase < np; ++phase) {
1651 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1655 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1656 const double dt = simulator.timeStepSize();
1658 bool converged =
false;
1659 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1660 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1662 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1666 well_potentials.clear();
1667 well_potentials.resize(np, 0.0);
1668 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1669 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1670 const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
1671 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1675 if constexpr (has_solvent) {
1677 assert(pu.phase_used[Gas]);
1678 const int gas_pos = pu.phase_pos[Gas];
1679 const EvalWell rate = well_copy.primary_variables_.getQs(Indices::contiSolventEqIdx);
1680 well_potentials[gas_pos] += rate.value();
1686 template<
typename TypeTag>
1687 typename StandardWell<TypeTag>::Scalar
1688 StandardWell<TypeTag>::
1689 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &simulator,
1690 const SummaryState &summary_state,
1691 DeferredLogger& deferred_logger,
1692 std::vector<Scalar>& potentials,
1696 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1697 simulator, summary_state, alq, deferred_logger,
true);
1698 if (bhp_at_thp_limit) {
1699 const auto& controls = this->well_ecl_.productionControls(summary_state);
1700 bhp = std::max(*bhp_at_thp_limit,
1701 static_cast<Scalar
>(controls.bhp_limit));
1702 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1705 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1706 "Failed in getting converged thp based potential calculation for well "
1707 + name() +
". Instead the bhp based value is used");
1708 const auto& controls = this->well_ecl_.productionControls(summary_state);
1709 bhp = controls.bhp_limit;
1710 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1715 template<
typename TypeTag>
1717 StandardWell<TypeTag>::
1718 computeWellRatesWithThpAlqProd(
const Simulator& simulator,
1719 const SummaryState& summary_state,
1720 DeferredLogger& deferred_logger,
1721 std::vector<Scalar>& potentials,
1725 computeWellRatesAndBhpWithThpAlqProd(simulator,
1732 template<
typename TypeTag>
1737 std::vector<Scalar>& well_potentials,
1740 const auto [compute_potential, bhp_controlled_well] =
1743 if (!compute_potential) {
1747 bool converged_implicit =
false;
1751 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1752 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1754 if (!converged_implicit) {
1756 const auto& summaryState = simulator.vanguard().summaryState();
1757 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1767 const auto& ws = well_state.well(this->index_of_well_);
1768 if (this->isInjector())
1769 bhp = std::max(ws.bhp, bhp);
1771 bhp = std::min(ws.bhp, bhp);
1773 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1774 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1777 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1781 this->checkNegativeWellPotentials(well_potentials,
1782 this->param_.check_well_operability_,
1792 template<
typename TypeTag>
1793 typename StandardWell<TypeTag>::Scalar
1796 const int openConnIdx)
const
1798 return (openConnIdx < 0)
1800 : this->connections_.rho(openConnIdx);
1807 template<
typename TypeTag>
1809 StandardWell<TypeTag>::
1810 updatePrimaryVariables(
const Simulator& simulator,
1811 const WellState<Scalar>& well_state,
1812 DeferredLogger& deferred_logger)
1814 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1816 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1817 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1820 if constexpr (Base::has_polymermw) {
1821 this->primary_variables_.updatePolyMW(well_state);
1824 this->primary_variables_.checkFinite(deferred_logger);
1830 template<
typename TypeTag>
1831 typename StandardWell<TypeTag>::Scalar
1832 StandardWell<TypeTag>::
1833 getRefDensity()
const
1835 return this->connections_.rho();
1841 template<
typename TypeTag>
1843 StandardWell<TypeTag>::
1844 updateWaterMobilityWithPolymer(
const Simulator& simulator,
1846 std::vector<EvalWell>& mob,
1847 DeferredLogger& deferred_logger)
const
1849 const int cell_idx = this->well_cells_[perf];
1850 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1851 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1855 if (this->isInjector()) {
1857 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1858 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1859 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1862 if (PolymerModule::hasPlyshlog()) {
1865 if (this->isInjector() && this->wpolymer() == 0.) {
1870 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1871 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1873 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1874 PerforationRates<Scalar> perf_rates;
1875 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1876 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1877 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1878 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1879 perf_rates, deferred_logger);
1881 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1882 const auto& material_law_manager = simulator.problem().materialLawManager();
1883 const auto& scaled_drainage_info =
1884 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1885 const Scalar swcr = scaled_drainage_info.Swcr;
1886 const EvalWell poro = this->extendEval(int_quant.porosity());
1887 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1889 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1890 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1891 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1893 if (PolymerModule::hasShrate()) {
1896 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1898 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1899 int_quant.pvtRegionIndex(),
1902 mob[waterCompIdx] /= shear_factor;
1906 template<
typename TypeTag>
1908 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
1910 this->linSys_.extract(jacobian);
1914 template <
typename TypeTag>
1916 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1917 const BVector& weights,
1918 const int pressureVarIndex,
1919 const bool use_well_weights,
1920 const WellState<Scalar>& well_state)
const
1922 this->linSys_.extractCPRPressureMatrix(jacobian,
1933 template<
typename TypeTag>
1934 typename StandardWell<TypeTag>::EvalWell
1935 StandardWell<TypeTag>::
1936 pskinwater(
const Scalar throughput,
1937 const EvalWell& water_velocity,
1938 DeferredLogger& deferred_logger)
const
1940 if constexpr (Base::has_polymermw) {
1941 const int water_table_id = this->polymerWaterTable_();
1942 if (water_table_id <= 0) {
1943 OPM_DEFLOG_THROW(std::runtime_error,
1944 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1947 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1948 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1950 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1951 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1954 OPM_DEFLOG_THROW(std::runtime_error,
1955 fmt::format(
"Polymermw is not activated, while injecting "
1956 "skin pressure is requested for well {}", name()),
1965 template<
typename TypeTag>
1966 typename StandardWell<TypeTag>::EvalWell
1967 StandardWell<TypeTag>::
1968 pskin(
const Scalar throughput,
1969 const EvalWell& water_velocity,
1970 const EvalWell& poly_inj_conc,
1971 DeferredLogger& deferred_logger)
const
1973 if constexpr (Base::has_polymermw) {
1974 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1975 const EvalWell water_velocity_abs = abs(water_velocity);
1976 if (poly_inj_conc == 0.) {
1977 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1979 const int polymer_table_id = this->polymerTable_();
1980 if (polymer_table_id <= 0) {
1981 OPM_DEFLOG_THROW(std::runtime_error,
1982 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1985 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1986 const Scalar reference_concentration = skprpolytable.refConcentration;
1987 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1989 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1990 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1991 if (poly_inj_conc == reference_concentration) {
1992 return sign * pskin_poly;
1995 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1996 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1997 return sign * pskin;
1999 OPM_DEFLOG_THROW(std::runtime_error,
2000 fmt::format(
"Polymermw is not activated, while injecting "
2001 "skin pressure is requested for well {}", name()),
2010 template<
typename TypeTag>
2011 typename StandardWell<TypeTag>::EvalWell
2012 StandardWell<TypeTag>::
2013 wpolymermw(
const Scalar throughput,
2014 const EvalWell& water_velocity,
2015 DeferredLogger& deferred_logger)
const
2017 if constexpr (Base::has_polymermw) {
2018 const int table_id = this->polymerInjTable_();
2019 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2020 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2021 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2022 if (this->wpolymer() == 0.) {
2023 return molecular_weight;
2025 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2026 return molecular_weight;
2028 OPM_DEFLOG_THROW(std::runtime_error,
2029 fmt::format(
"Polymermw is not activated, while injecting "
2030 "polymer molecular weight is requested for well {}", name()),
2039 template<
typename TypeTag>
2041 StandardWell<TypeTag>::
2042 updateWaterThroughput([[maybe_unused]]
const double dt,
2043 WellState<Scalar>& well_state)
const
2045 if constexpr (Base::has_polymermw) {
2046 if (!this->isInjector()) {
2050 auto& perf_water_throughput = well_state.well(this->index_of_well_)
2051 .perf_data.water_throughput;
2053 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2054 const Scalar perf_water_vel =
2055 this->primary_variables_.value(Bhp + 1 + perf);
2059 if (perf_water_vel > Scalar{0}) {
2060 perf_water_throughput[perf] += perf_water_vel * dt;
2070 template<
typename TypeTag>
2072 StandardWell<TypeTag>::
2073 handleInjectivityRate(
const Simulator& simulator,
2075 std::vector<EvalWell>& cq_s)
const
2077 const int cell_idx = this->well_cells_[perf];
2078 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2079 const auto& fs = int_quants.fluidState();
2080 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2081 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2082 const int wat_vel_index = Bhp + 1 + perf;
2083 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2087 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2093 template<
typename TypeTag>
2095 StandardWell<TypeTag>::
2096 handleInjectivityEquations(
const Simulator& simulator,
2097 const WellState<Scalar>& well_state,
2099 const EvalWell& water_flux_s,
2100 DeferredLogger& deferred_logger)
2102 const int cell_idx = this->well_cells_[perf];
2103 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2104 const auto& fs = int_quants.fluidState();
2105 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2106 const EvalWell water_flux_r = water_flux_s / b_w;
2107 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2108 const EvalWell water_velocity = water_flux_r / area;
2109 const int wat_vel_index = Bhp + 1 + perf;
2112 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2114 const auto& ws = well_state.well(this->index_of_well_);
2115 const auto& perf_data = ws.perf_data;
2116 const auto& perf_water_throughput = perf_data.water_throughput;
2117 const Scalar throughput = perf_water_throughput[perf];
2118 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2120 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2121 poly_conc.setValue(this->wpolymer());
2124 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2125 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2127 StandardWellAssemble<FluidSystem,Indices>(*this).
2128 assembleInjectivityEq(eq_pskin,
2133 this->primary_variables_.numWellEq(),
2141 template<
typename TypeTag>
2143 StandardWell<TypeTag>::
2144 checkConvergenceExtraEqs(
const std::vector<Scalar>& res,
2145 ConvergenceReport& report)
const
2150 if constexpr (Base::has_polymermw) {
2151 WellConvergence(*this).
2152 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2160 template<
typename TypeTag>
2162 StandardWell<TypeTag>::
2163 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2164 const IntensiveQuantities& int_quants,
2165 const WellState<Scalar>& well_state,
2167 std::vector<RateVector>& connectionRates,
2168 DeferredLogger& deferred_logger)
const
2171 EvalWell cq_s_polymw = cq_s_poly;
2172 if (this->isInjector()) {
2173 const int wat_vel_index = Bhp + 1 + perf;
2174 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2175 if (water_velocity > 0.) {
2176 const auto& ws = well_state.well(this->index_of_well_);
2177 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2178 const Scalar throughput = perf_water_throughput[perf];
2179 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2180 cq_s_polymw *= molecular_weight;
2186 }
else if (this->isProducer()) {
2187 if (cq_s_polymw < 0.) {
2188 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2195 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2203 template<
typename TypeTag>
2204 std::optional<typename StandardWell<TypeTag>::Scalar>
2205 StandardWell<TypeTag>::
2206 computeBhpAtThpLimitProd(
const WellState<Scalar>& well_state,
2207 const Simulator& simulator,
2208 const SummaryState& summary_state,
2209 DeferredLogger& deferred_logger)
const
2211 return computeBhpAtThpLimitProdWithAlq(simulator,
2213 this->getALQ(well_state),
2218 template<
typename TypeTag>
2219 std::optional<typename StandardWell<TypeTag>::Scalar>
2220 StandardWell<TypeTag>::
2221 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2222 const SummaryState& summary_state,
2223 const Scalar alq_value,
2224 DeferredLogger& deferred_logger,
2225 bool iterate_if_no_solution)
const
2229 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2235 std::vector<Scalar> rates(3);
2236 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2237 this->adaptRatesForVFP(rates);
2241 Scalar max_pressure = 0.0;
2242 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2243 const int cell_idx = this->well_cells_[perf];
2244 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2245 const auto& fs = int_quants.fluidState();
2246 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2247 max_pressure = std::max(max_pressure, pressure_cell);
2249 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2252 this->connections_.rho(),
2254 this->getTHPConstraint(summary_state),
2258 auto v = frates(*bhpAtLimit);
2259 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2264 if (!iterate_if_no_solution)
2265 return std::nullopt;
2267 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2271 std::vector<Scalar> rates(3);
2272 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2273 this->adaptRatesForVFP(rates);
2277 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2280 this->connections_.rho(),
2282 this->getTHPConstraint(summary_state),
2288 auto v = frates(*bhpAtLimit);
2289 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2295 return std::nullopt;
2300 template<
typename TypeTag>
2301 std::optional<typename StandardWell<TypeTag>::Scalar>
2302 StandardWell<TypeTag>::
2303 computeBhpAtThpLimitInj(
const Simulator& simulator,
2304 const SummaryState& summary_state,
2305 DeferredLogger& deferred_logger)
const
2308 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2314 std::vector<Scalar> rates(3);
2315 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2319 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2321 this->connections_.rho(),
2332 template<
typename TypeTag>
2334 StandardWell<TypeTag>::
2335 iterateWellEqWithControl(
const Simulator& simulator,
2337 const Well::InjectionControls& inj_controls,
2338 const Well::ProductionControls& prod_controls,
2339 WellState<Scalar>& well_state,
2340 const GroupState<Scalar>& group_state,
2341 DeferredLogger& deferred_logger)
2343 const int max_iter = this->param_.max_inner_iter_wells_;
2346 bool relax_convergence =
false;
2347 this->regularize_ =
false;
2349 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2351 if (it > this->param_.strict_inner_iter_wells_) {
2352 relax_convergence =
true;
2353 this->regularize_ =
true;
2356 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2358 converged = report.converged();
2364 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2371 }
while (it < max_iter);
2377 template<
typename TypeTag>
2379 StandardWell<TypeTag>::
2380 iterateWellEqWithSwitching(
const Simulator& simulator,
2382 const Well::InjectionControls& inj_controls,
2383 const Well::ProductionControls& prod_controls,
2384 WellState<Scalar>& well_state,
2385 const GroupState<Scalar>& group_state,
2386 DeferredLogger& deferred_logger,
2387 const bool fixed_control ,
2388 const bool fixed_status )
2390 const int max_iter = this->param_.max_inner_iter_wells_;
2392 bool converged =
false;
2393 bool relax_convergence =
false;
2394 this->regularize_ =
false;
2395 const auto& summary_state = simulator.vanguard().summaryState();
2400 constexpr int min_its_after_switch = 4;
2401 int its_since_last_switch = min_its_after_switch;
2402 int switch_count= 0;
2404 const auto well_status_orig = this->wellStatus_;
2405 const auto operability_orig = this->operability_status_;
2406 auto well_status_cur = well_status_orig;
2407 int status_switch_count = 0;
2409 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2410 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2412 const bool allow_switching =
2413 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2414 (!fixed_control || !fixed_status) && allow_open;
2416 bool changed =
false;
2417 bool final_check =
false;
2419 this->operability_status_.resetOperability();
2420 this->operability_status_.solvable =
true;
2422 its_since_last_switch++;
2423 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2424 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2425 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2426 inj_controls, prod_controls, wqTotal,
2427 deferred_logger, fixed_control, fixed_status);
2429 its_since_last_switch = 0;
2431 if (well_status_cur != this->wellStatus_) {
2432 well_status_cur = this->wellStatus_;
2433 status_switch_count++;
2436 if (!changed && final_check) {
2439 final_check =
false;
2443 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2445 if (it > this->param_.strict_inner_iter_wells_) {
2446 relax_convergence =
true;
2447 this->regularize_ =
true;
2450 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2452 converged = report.converged();
2456 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2458 its_since_last_switch = min_its_after_switch;
2465 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2467 }
while (it < max_iter);
2470 if (allow_switching){
2472 const bool is_stopped = this->wellIsStopped();
2473 if (this->wellHasTHPConstraints(summary_state)){
2474 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2475 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2477 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2481 this->wellStatus_ = well_status_orig;
2482 this->operability_status_ = operability_orig;
2483 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2484 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2485 deferred_logger.debug(message);
2491 template<
typename TypeTag>
2492 std::vector<typename StandardWell<TypeTag>::Scalar>
2498 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2499 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2500 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2501 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2502 const int cell_idx = this->well_cells_[perf];
2503 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2504 std::vector<Scalar> mob(this->num_components_, 0.);
2505 getMobility(simulator, perf, mob, deferred_logger);
2506 std::vector<Scalar> cq_s(this->num_components_, 0.);
2507 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2508 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2509 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2511 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2512 cq_s, perf_rates, deferred_logger);
2513 for (
int comp = 0; comp < this->num_components_; ++comp) {
2514 well_q_s[comp] += cq_s[comp];
2517 const auto& comm = this->parallel_well_info_.communication();
2518 if (comm.size() > 1)
2520 comm.sum(well_q_s.data(), well_q_s.size());
2527 template <
typename TypeTag>
2528 std::vector<typename StandardWell<TypeTag>::Scalar>
2532 const int num_pri_vars = this->primary_variables_.numWellEq();
2533 std::vector<Scalar> retval(num_pri_vars);
2534 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2535 retval[ii] = this->primary_variables_.value(ii);
2544 template <
typename TypeTag>
2546 StandardWell<TypeTag>::
2547 setPrimaryVars(
typename std::vector<Scalar>::const_iterator it)
2549 const int num_pri_vars = this->primary_variables_.numWellEq();
2550 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2551 this->primary_variables_.setValue(ii, it[ii]);
2553 return num_pri_vars;
2557 template <
typename TypeTag>
2558 typename StandardWell<TypeTag>::Eval
2559 StandardWell<TypeTag>::
2560 connectionRateEnergy(
const Scalar maxOilSaturation,
2561 const std::vector<EvalWell>& cq_s,
2562 const IntensiveQuantities& intQuants,
2563 DeferredLogger& deferred_logger)
const
2565 auto fs = intQuants.fluidState();
2567 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2568 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2573 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2574 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2575 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2576 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2577 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2580 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2581 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2586 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2588 deferred_logger.debug(
2589 fmt::format(
"Problematic d value {} obtained for well {}"
2590 " during calculateSinglePerf with rs {}"
2591 ", rv {}. Continue as if no dissolution (rs = 0) and"
2592 " vaporization (rv = 0) for this connection.",
2593 d, this->name(), fs.Rs(), fs.Rv()));
2594 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2596 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2597 cq_r_thermal = (cq_s[gasCompIdx] -
2598 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2599 (d * this->extendEval(fs.invB(phaseIdx)) );
2600 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2602 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2604 (d * this->extendEval(fs.invB(phaseIdx)) );
2610 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2612 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2613 fs.setTemperature(this->well_ecl_.inj_temperature());
2614 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2615 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2616 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2617 paramCache.setRegionIndex(pvtRegionIdx);
2618 paramCache.setMaxOilSat(maxOilSaturation);
2619 paramCache.updatePhase(fs, phaseIdx);
2621 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2622 fs.setDensity(phaseIdx, rho);
2623 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2624 fs.setEnthalpy(phaseIdx, h);
2625 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2626 result += getValue(cq_r_thermal);
2627 }
else if (cq_r_thermal > 0.0) {
2628 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2629 result += Base::restrictEval(cq_r_thermal);
2632 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2633 result += Base::restrictEval(cq_r_thermal);
2637 return result * this->well_efficiency_factor_;
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
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:304
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:291
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:1413
Class for computing BHP limits.
Definition WellBhpThpCalculator.hpp:41
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
Definition WellBhpThpCalculator.cpp:93
Definition WellInterfaceGeneric.hpp:53
Collect per-connection static information to enable calculating connection-level or well-level produc...
Definition WellProdIndexCalculator.hpp:37
Scalar connectionProdIndStandard(const std::size_t connIdx, const Scalar connMobility) const
Compute connection-level steady-state productivity index value using dynamic phase mobility.
Definition WellProdIndexCalculator.cpp:121
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:66
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition PerforationData.hpp:41