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(perf_data.microbial_rates[perf],
627 perf_data.oxygen_rates[perf],
628 perf_data.urea_rates[perf],
630 microbialConcentration,
636 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
639 const auto& pu = well_state.phaseUsage();
640 if (pu.has_co2_or_h2store) {
641 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
642 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
643 perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
649 template<
typename TypeTag>
650 template<
class Value>
652 StandardWell<TypeTag>::
653 getMobility(
const Simulator& simulator,
655 std::vector<Value>& mob,
656 DeferredLogger& deferred_logger)
const
658 auto obtain = [
this](
const Eval& value)
660 if constexpr (std::is_same_v<Value, Scalar>) {
661 static_cast<void>(
this);
662 return getValue(value);
664 return this->extendEval(value);
667 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
668 obtain, deferred_logger);
671 if constexpr (has_polymer) {
672 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
673 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
678 if constexpr (!Base::has_polymermw) {
679 if constexpr (std::is_same_v<Value, Scalar>) {
680 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
681 for (std::size_t i = 0; i < mob.size(); ++i) {
682 mob_eval[i].setValue(mob[i]);
684 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
685 for (std::size_t i = 0; i < mob.size(); ++i) {
686 mob[i] = getValue(mob_eval[i]);
689 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
695 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
696 const Scalar bhp = this->primary_variables_.value(Bhp);
697 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
698 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
699 for (std::size_t i = 0; i < mob.size(); ++i) {
700 mob[i] *= multiplier;
706 template<
typename TypeTag>
708 StandardWell<TypeTag>::
709 updateWellState(
const Simulator& simulator,
710 const BVectorWell& dwells,
711 WellState<Scalar>& well_state,
712 DeferredLogger& deferred_logger)
714 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
716 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
717 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
719 const auto& summary_state = simulator.vanguard().summaryState();
720 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
721 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
728 template<
typename TypeTag>
730 StandardWell<TypeTag>::
731 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
732 const bool stop_or_zero_rate_target,
733 DeferredLogger& deferred_logger)
735 const Scalar dFLimit = this->param_.dwell_fraction_max_;
736 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
737 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
740 if constexpr (Base::has_polymermw) {
741 this->primary_variables_.updateNewtonPolyMW(dwells);
744 this->primary_variables_.checkFinite(deferred_logger);
751 template<
typename TypeTag>
753 StandardWell<TypeTag>::
754 updateWellStateFromPrimaryVariables(WellState<Scalar>& well_state,
755 const SummaryState& summary_state,
756 DeferredLogger& deferred_logger)
const
758 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
761 if constexpr (Base::has_polymermw) {
762 this->primary_variables_.copyToWellStatePolyMW(well_state);
770 template<
typename TypeTag>
772 StandardWell<TypeTag>::
773 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
778 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
779 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
781 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
782 std::vector<Scalar> mob(this->num_components_, 0.0);
783 getMobility(simulator, perf, mob, deferred_logger);
785 const int cell_idx = this->well_cells_[perf];
786 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
787 const auto& fs = int_quantities.fluidState();
789 Scalar p_r = this->getPerfCellPressure(fs).value();
792 std::vector<Scalar> b_perf(this->num_components_);
793 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
794 if (!FluidSystem::phaseIsActive(phase)) {
797 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
798 b_perf[comp_idx] = fs.invB(phase).value();
800 if constexpr (has_solvent) {
801 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
805 const Scalar h_perf = this->connections_.pressure_diff(perf);
806 const Scalar pressure_diff = p_r - h_perf;
811 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
812 deferred_logger.debug(
"CROSSFLOW_IPR",
813 "cross flow found when updateIPR for well " + name()
814 +
" . The connection is ignored in IPR calculations");
820 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
821 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
822 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
826 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
827 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
828 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
829 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
830 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
831 ipr_b_perf[comp_idx] += tw_mob;
835 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
836 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
837 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
838 const Scalar rs = (fs.Rs()).value();
839 const Scalar rv = (fs.Rv()).value();
841 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
842 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
844 ipr_a_perf[gas_comp_idx] += dis_gas_a;
845 ipr_a_perf[oil_comp_idx] += vap_oil_a;
847 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
848 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
850 ipr_b_perf[gas_comp_idx] += dis_gas_b;
851 ipr_b_perf[oil_comp_idx] += vap_oil_b;
854 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
855 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
856 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
859 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
860 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
863 template<
typename TypeTag>
865 StandardWell<TypeTag>::
866 updateIPRImplicit(
const Simulator& simulator,
867 WellState<Scalar>& well_state,
868 DeferredLogger& deferred_logger)
877 auto rates = well_state.well(this->index_of_well_).surface_rates;
879 for (std::size_t p = 0; p < rates.size(); ++p) {
880 zero_rates &= rates[p] == 0.0;
882 auto& ws = well_state.well(this->index_of_well_);
884 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
885 deferred_logger.debug(msg);
897 const auto& group_state = simulator.problem().wellModel().groupState();
899 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
900 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
902 auto inj_controls = Well::InjectionControls(0);
903 auto prod_controls = Well::ProductionControls(0);
904 prod_controls.addControl(Well::ProducerCMode::BHP);
905 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
908 const auto cmode = ws.production_cmode;
909 ws.production_cmode = Well::ProducerCMode::BHP;
910 const double dt = simulator.timeStepSize();
911 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
913 const size_t nEq = this->primary_variables_.numWellEq();
917 for (
size_t i=0; i < nEq; ++i){
922 BVectorWell x_well(1);
923 x_well[0].resize(nEq);
924 this->linSys_.solve(rhs, x_well);
926 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
927 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
928 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
929 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
931 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
933 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
936 ws.production_cmode = cmode;
939 template<
typename TypeTag>
941 StandardWell<TypeTag>::
942 checkOperabilityUnderBHPLimit(
const WellState<Scalar>& well_state,
943 const Simulator& simulator,
944 DeferredLogger& deferred_logger)
946 const auto& summaryState = simulator.vanguard().summaryState();
947 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
950 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
951 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
954 Scalar total_ipr_mass_rate = 0.0;
955 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
957 if (!FluidSystem::phaseIsActive(phaseIdx)) {
961 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
962 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
964 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
965 total_ipr_mass_rate += ipr_rate * rho;
967 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
968 this->operability_status_.operable_under_only_bhp_limit =
false;
972 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
976 std::vector<Scalar> well_rates_bhp_limit;
977 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
979 this->adaptRatesForVFP(well_rates_bhp_limit);
980 const Scalar thp_limit = this->getTHPConstraint(summaryState);
981 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
983 this->connections_.rho(),
984 this->getALQ(well_state),
987 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
988 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
999 this->operability_status_.operable_under_only_bhp_limit =
true;
1000 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1008 template<
typename TypeTag>
1010 StandardWell<TypeTag>::
1011 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
1012 const WellState<Scalar>& well_state,
1013 DeferredLogger& deferred_logger)
1015 const auto& summaryState = simulator.vanguard().summaryState();
1016 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1017 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1020 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1022 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1023 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1024 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1026 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1027 if (this->isProducer() && *obtain_bhp < thp_limit) {
1028 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1029 +
" bars is SMALLER than thp limit "
1030 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1031 +
" bars as a producer for well " + name();
1032 deferred_logger.debug(msg);
1034 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1035 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1036 +
" bars is LARGER than thp limit "
1037 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1038 +
" bars as a injector for well " + name();
1039 deferred_logger.debug(msg);
1042 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1043 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1044 if (!this->wellIsStopped()) {
1045 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1046 deferred_logger.debug(
" could not find bhp value at thp limit "
1047 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1048 +
" bar for well " + name() +
", the well might need to be closed ");
1057 template<
typename TypeTag>
1059 StandardWell<TypeTag>::
1060 allDrawDownWrongDirection(
const Simulator& simulator)
const
1062 bool all_drawdown_wrong_direction =
true;
1064 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1065 const int cell_idx = this->well_cells_[perf];
1066 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1067 const auto& fs = intQuants.fluidState();
1069 const Scalar pressure = this->getPerfCellPressure(fs).value();
1070 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1073 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1074 const Scalar drawdown = pressure - well_pressure;
1079 if ( (drawdown < 0. && this->isInjector()) ||
1080 (drawdown > 0. && this->isProducer()) ) {
1081 all_drawdown_wrong_direction =
false;
1086 const auto& comm = this->parallel_well_info_.communication();
1087 if (comm.size() > 1)
1089 all_drawdown_wrong_direction =
1090 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1093 return all_drawdown_wrong_direction;
1099 template<
typename TypeTag>
1101 StandardWell<TypeTag>::
1102 canProduceInjectWithCurrentBhp(
const Simulator& simulator,
1103 const WellState<Scalar>& well_state,
1104 DeferredLogger& deferred_logger)
1106 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
1107 std::vector<Scalar> well_rates;
1108 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1110 const Scalar sign = (this->isProducer()) ? -1. : 1.;
1111 const Scalar threshold = sign * std::numeric_limits<Scalar>::min();
1113 bool can_produce_inject =
false;
1114 for (
const auto value : well_rates) {
1115 if (this->isProducer() && value < threshold) {
1116 can_produce_inject =
true;
1118 }
else if (this->isInjector() && value > threshold) {
1119 can_produce_inject =
true;
1124 if (!can_produce_inject) {
1125 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1128 return can_produce_inject;
1135 template<
typename TypeTag>
1137 StandardWell<TypeTag>::
1138 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1140 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1146 template<
typename TypeTag>
1147 typename StandardWell<TypeTag>::WellConnectionProps
1148 StandardWell<TypeTag>::
1149 computePropertiesForWellConnectionPressures(
const Simulator& simulator,
1150 const WellState<Scalar>& well_state)
const
1152 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1154 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1156 return model.intensiveQuantities(cell_idx, 0)
1157 .fluidState().temperature(phase_idx).value();
1161 [&model = simulator.model()](
int cell_idx)
1163 return model.intensiveQuantities(cell_idx, 0)
1164 .fluidState().saltConcentration().value();
1168 [&model = simulator.model()](
int cell_idx)
1170 return model.intensiveQuantities(cell_idx, 0)
1171 .fluidState().pvtRegionIndex();
1175 if constexpr (Indices::enableSolvent) {
1176 prop_func.solventInverseFormationVolumeFactor =
1177 [&model = simulator.model()](
int cell_idx)
1179 return model.intensiveQuantities(cell_idx, 0)
1180 .solventInverseFormationVolumeFactor().value();
1183 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1185 return model.intensiveQuantities(cell_idx, 0)
1186 .solventRefDensity();
1190 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1197 template<
typename TypeTag>
1202 const std::vector<Scalar>& B_avg,
1204 const bool relax_tolerance)
const
1208 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1210 Scalar tol_wells = this->param_.tolerance_wells_;
1212 constexpr Scalar stopped_factor = 1.e-4;
1214 constexpr Scalar dynamic_thp_factor = 1.e-1;
1215 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1216 tol_wells = tol_wells*stopped_factor;
1217 }
else if (this->getDynamicThpLimit()) {
1218 tol_wells = tol_wells*dynamic_thp_factor;
1221 std::vector<Scalar> res;
1224 this->param_.max_residual_allowed_,
1226 this->param_.relaxed_tolerance_flow_well_,
1228 this->wellIsStopped(),
1232 checkConvergenceExtraEqs(res, report);
1241 template<
typename TypeTag>
1249 auto fluidState = [&simulator,
this](
const int perf)
1251 const auto cell_idx = this->well_cells_[perf];
1252 return simulator.
model()
1253 .intensiveQuantities(cell_idx, 0).fluidState();
1256 const int np = this->number_of_phases_;
1257 auto setToZero = [np](Scalar* x) ->
void
1259 std::fill_n(x, np, 0.0);
1262 auto addVector = [np](
const Scalar* src, Scalar* dest) ->
void
1264 std::transform(src, src + np, dest, dest, std::plus<>{});
1267 auto& ws = well_state.well(this->index_of_well_);
1268 auto& perf_data = ws.perf_data;
1269 auto* wellPI = ws.productivity_index.data();
1270 auto* connPI = perf_data.prod_index.data();
1274 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1275 auto subsetPerfID = 0;
1277 for (
const auto& perf : *this->perf_data_) {
1278 auto allPerfID = perf.ecl_index;
1280 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) -> Scalar
1285 std::vector<Scalar> mob(this->num_components_, 0.0);
1286 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1288 const auto& fs = fluidState(subsetPerfID);
1291 if (this->isInjector()) {
1292 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1293 mob, connPI, deferred_logger);
1296 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1299 addVector(connPI, wellPI);
1306 const auto& comm = this->parallel_well_info_.communication();
1307 if (comm.size() > 1) {
1308 comm.sum(wellPI, np);
1311 assert ((
static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1312 "Internal logic error in processing connections for PI/II");
1317 template<
typename TypeTag>
1318 void StandardWell<TypeTag>::
1319 computeWellConnectionDensitesPressures(
const Simulator& simulator,
1320 const WellState<Scalar>& well_state,
1321 const WellConnectionProps& props,
1322 DeferredLogger& deferred_logger)
1327 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1332 [&model = simulator.model()](
const int cell,
1333 const std::vector<int>& phases,
1334 std::vector<Scalar>& mob)
1336 const auto& iq = model.intensiveQuantities(cell, 0);
1338 std::transform(phases.begin(), phases.end(), mob.begin(),
1339 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1344 [&model = simulator.model()](
const int cell,
1345 const std::vector<int>& phases,
1346 std::vector<Scalar>& rho)
1348 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1350 std::transform(phases.begin(), phases.end(), rho.begin(),
1351 [&fs](
const int phase) { return fs.density(phase).value(); });
1355 const auto stopped_or_zero_rate_target = this->
1356 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1359 .computeProperties(stopped_or_zero_rate_target, well_state,
1360 prop_func, props, deferred_logger);
1367 template<
typename TypeTag>
1369 StandardWell<TypeTag>::
1370 computeWellConnectionPressures(
const Simulator& simulator,
1371 const WellState<Scalar>& well_state,
1372 DeferredLogger& deferred_logger)
1374 const auto props = computePropertiesForWellConnectionPressures
1375 (simulator, well_state);
1377 computeWellConnectionDensitesPressures(simulator, well_state,
1378 props, deferred_logger);
1385 template<
typename TypeTag>
1387 StandardWell<TypeTag>::
1388 solveEqAndUpdateWellState(
const Simulator& simulator,
1389 WellState<Scalar>& well_state,
1390 DeferredLogger& deferred_logger)
1392 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1396 BVectorWell dx_well(1);
1397 dx_well[0].resize(this->primary_variables_.numWellEq());
1398 this->linSys_.solve( dx_well);
1400 updateWellState(simulator, dx_well, well_state, deferred_logger);
1407 template<
typename TypeTag>
1409 StandardWell<TypeTag>::
1410 calculateExplicitQuantities(
const Simulator& simulator,
1411 const WellState<Scalar>& well_state,
1412 DeferredLogger& deferred_logger)
1414 updatePrimaryVariables(simulator, well_state, deferred_logger);
1415 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1416 this->computeAccumWell();
1421 template<
typename TypeTag>
1424 apply(
const BVector& x, BVector& Ax)
const
1426 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1428 if (this->param_.matrix_add_well_contributions_)
1434 this->linSys_.
apply(x, Ax);
1440 template<
typename TypeTag>
1443 apply(BVector& r)
const
1445 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1447 this->linSys_.
apply(r);
1453 template<
typename TypeTag>
1461 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1464 xw[0].resize(this->primary_variables_.numWellEq());
1466 this->linSys_.recoverSolutionWell(x, xw);
1467 updateWellState(simulator, xw, well_state, deferred_logger);
1473 template<
typename TypeTag>
1478 std::vector<Scalar>& well_flux,
1482 const int np = this->number_of_phases_;
1483 well_flux.resize(np, 0.0);
1485 const bool allow_cf = this->getAllowCrossFlow();
1487 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1488 const int cell_idx = this->well_cells_[perf];
1489 const auto& intQuants = simulator.
model().intensiveQuantities(cell_idx, 0);
1491 std::vector<Scalar> mob(this->num_components_, 0.);
1492 getMobility(simulator, perf, mob, deferred_logger);
1493 Scalar trans_mult = simulator.
problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1494 const auto& wellstate_nupcol = simulator.
problem().wellModel().nupcolWellState().well(this->index_of_well_);
1495 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1497 std::vector<Scalar> cq_s(this->num_components_, 0.);
1499 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1500 cq_s, perf_rates, deferred_logger);
1502 for(
int p = 0; p < np; ++p) {
1503 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1507 if constexpr (has_solvent) {
1508 const auto& pu = this->phaseUsage();
1509 assert(pu.phase_used[Gas]);
1510 const int gas_pos = pu.phase_pos[Gas];
1511 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1514 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1519 template<
typename TypeTag>
1521 StandardWell<TypeTag>::
1522 computeWellRatesWithBhpIterations(
const Simulator& simulator,
1524 std::vector<Scalar>& well_flux,
1525 DeferredLogger& deferred_logger)
const
1529 StandardWell<TypeTag> well_copy(*
this);
1530 well_copy.resetDampening();
1535 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1536 const auto& group_state = simulator.problem().wellModel().groupState();
1539 const auto& summary_state = simulator.vanguard().summaryState();
1540 auto inj_controls = well_copy.well_ecl_.isInjector()
1541 ? well_copy.well_ecl_.injectionControls(summary_state)
1542 : Well::InjectionControls(0);
1543 auto prod_controls = well_copy.well_ecl_.isProducer()
1544 ? well_copy.well_ecl_.productionControls(summary_state) :
1545 Well::ProductionControls(0);
1548 auto& ws = well_state_copy.well(this->index_of_well_);
1549 if (well_copy.well_ecl_.isInjector()) {
1550 inj_controls.bhp_limit = bhp;
1551 ws.injection_cmode = Well::InjectorCMode::BHP;
1553 prod_controls.bhp_limit = bhp;
1554 ws.production_cmode = Well::ProducerCMode::BHP;
1559 const int np = this->number_of_phases_;
1560 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1561 for (
int phase = 0; phase < np; ++phase){
1562 well_state_copy.wellRates(this->index_of_well_)[phase]
1563 = sign * ws.well_potentials[phase];
1565 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1566 well_copy.computeAccumWell();
1568 const double dt = simulator.timeStepSize();
1569 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1571 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1572 " potentials are computed based on unconverged solution";
1573 deferred_logger.debug(msg);
1575 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1576 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1577 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1583 template<
typename TypeTag>
1584 std::vector<typename StandardWell<TypeTag>::Scalar>
1585 StandardWell<TypeTag>::
1586 computeWellPotentialWithTHP(
const Simulator& simulator,
1587 DeferredLogger& deferred_logger,
1588 const WellState<Scalar>& well_state)
const
1590 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1591 const auto& summary_state = simulator.vanguard().summaryState();
1593 const auto& well = this->well_ecl_;
1594 if (well.isInjector()){
1595 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1596 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1597 if (bhp_at_thp_limit) {
1598 const Scalar bhp = std::min(*bhp_at_thp_limit,
1599 static_cast<Scalar
>(controls.bhp_limit));
1600 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1602 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1603 "Failed in getting converged thp based potential calculation for well "
1604 + name() +
". Instead the bhp based value is used");
1605 const Scalar bhp = controls.bhp_limit;
1606 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1609 computeWellRatesWithThpAlqProd(
1610 simulator, summary_state,
1611 deferred_logger, potentials, this->getALQ(well_state)
1618 template<
typename TypeTag>
1620 StandardWell<TypeTag>::
1621 computeWellPotentialsImplicit(
const Simulator& simulator,
1622 const WellState<Scalar>& well_state,
1623 std::vector<Scalar>& well_potentials,
1624 DeferredLogger& deferred_logger)
const
1629 StandardWell<TypeTag> well_copy(*
this);
1632 WellState<Scalar> well_state_copy = well_state;
1633 const auto& group_state = simulator.problem().wellModel().groupState();
1634 auto& ws = well_state_copy.well(this->index_of_well_);
1637 const auto& summary_state = simulator.vanguard().summaryState();
1638 auto inj_controls = well_copy.well_ecl_.isInjector()
1639 ? well_copy.well_ecl_.injectionControls(summary_state)
1640 : Well::InjectionControls(0);
1641 auto prod_controls = well_copy.well_ecl_.isProducer()
1642 ? well_copy.well_ecl_.productionControls(summary_state) :
1643 Well::ProductionControls(0);
1646 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1649 const int num_perf = ws.perf_data.size();
1650 for (
int perf = 0; perf < num_perf; ++perf) {
1651 ws.perf_data.pressure[perf] = ws.bhp + well_copy.connections_.pressure_diff(perf);
1654 const int np = this->number_of_phases_;
1655 bool trivial =
true;
1656 for (
int phase = 0; phase < np; ++phase){
1657 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1660 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1661 for (
int phase = 0; phase < np; ++phase) {
1662 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1666 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1667 const double dt = simulator.timeStepSize();
1669 bool converged =
false;
1670 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1671 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1673 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1677 well_potentials.clear();
1678 well_potentials.resize(np, 0.0);
1679 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1680 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1681 const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
1682 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1686 if constexpr (has_solvent) {
1688 assert(pu.phase_used[Gas]);
1689 const int gas_pos = pu.phase_pos[Gas];
1690 const EvalWell rate = well_copy.primary_variables_.getQs(Indices::contiSolventEqIdx);
1691 well_potentials[gas_pos] += rate.value();
1697 template<
typename TypeTag>
1698 typename StandardWell<TypeTag>::Scalar
1699 StandardWell<TypeTag>::
1700 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &simulator,
1701 const SummaryState &summary_state,
1702 DeferredLogger& deferred_logger,
1703 std::vector<Scalar>& potentials,
1707 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1708 simulator, summary_state, alq, deferred_logger,
true);
1709 if (bhp_at_thp_limit) {
1710 const auto& controls = this->well_ecl_.productionControls(summary_state);
1711 bhp = std::max(*bhp_at_thp_limit,
1712 static_cast<Scalar
>(controls.bhp_limit));
1713 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1716 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1717 "Failed in getting converged thp based potential calculation for well "
1718 + name() +
". Instead the bhp based value is used");
1719 const auto& controls = this->well_ecl_.productionControls(summary_state);
1720 bhp = controls.bhp_limit;
1721 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1726 template<
typename TypeTag>
1728 StandardWell<TypeTag>::
1729 computeWellRatesWithThpAlqProd(
const Simulator& simulator,
1730 const SummaryState& summary_state,
1731 DeferredLogger& deferred_logger,
1732 std::vector<Scalar>& potentials,
1736 computeWellRatesAndBhpWithThpAlqProd(simulator,
1743 template<
typename TypeTag>
1748 std::vector<Scalar>& well_potentials,
1751 const auto [compute_potential, bhp_controlled_well] =
1754 if (!compute_potential) {
1758 bool converged_implicit =
false;
1762 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1763 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1765 if (!converged_implicit) {
1767 const auto& summaryState = simulator.vanguard().summaryState();
1768 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1778 const auto& ws = well_state.well(this->index_of_well_);
1779 if (this->isInjector())
1780 bhp = std::max(ws.bhp, bhp);
1782 bhp = std::min(ws.bhp, bhp);
1784 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1785 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1788 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1792 this->checkNegativeWellPotentials(well_potentials,
1793 this->param_.check_well_operability_,
1803 template<
typename TypeTag>
1804 typename StandardWell<TypeTag>::Scalar
1807 const int openConnIdx)
const
1809 return (openConnIdx < 0)
1811 : this->connections_.rho(openConnIdx);
1818 template<
typename TypeTag>
1820 StandardWell<TypeTag>::
1821 updatePrimaryVariables(
const Simulator& simulator,
1822 const WellState<Scalar>& well_state,
1823 DeferredLogger& deferred_logger)
1825 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1827 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1828 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1831 if constexpr (Base::has_polymermw) {
1832 this->primary_variables_.updatePolyMW(well_state);
1835 this->primary_variables_.checkFinite(deferred_logger);
1841 template<
typename TypeTag>
1842 typename StandardWell<TypeTag>::Scalar
1843 StandardWell<TypeTag>::
1844 getRefDensity()
const
1846 return this->connections_.rho();
1852 template<
typename TypeTag>
1854 StandardWell<TypeTag>::
1855 updateWaterMobilityWithPolymer(
const Simulator& simulator,
1857 std::vector<EvalWell>& mob,
1858 DeferredLogger& deferred_logger)
const
1860 const int cell_idx = this->well_cells_[perf];
1861 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1862 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1866 if (this->isInjector()) {
1868 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1869 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1870 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1873 if (PolymerModule::hasPlyshlog()) {
1876 if (this->isInjector() && this->wpolymer() == 0.) {
1881 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1882 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1884 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1885 PerforationRates<Scalar> perf_rates;
1886 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1887 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1888 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1889 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1890 perf_rates, deferred_logger);
1892 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1893 const auto& material_law_manager = simulator.problem().materialLawManager();
1894 const auto& scaled_drainage_info =
1895 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1896 const Scalar swcr = scaled_drainage_info.Swcr;
1897 const EvalWell poro = this->extendEval(int_quant.porosity());
1898 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1900 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1901 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1902 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1904 if (PolymerModule::hasShrate()) {
1907 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1909 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1910 int_quant.pvtRegionIndex(),
1913 mob[waterCompIdx] /= shear_factor;
1917 template<
typename TypeTag>
1919 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
1921 this->linSys_.extract(jacobian);
1925 template <
typename TypeTag>
1927 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1928 const BVector& weights,
1929 const int pressureVarIndex,
1930 const bool use_well_weights,
1931 const WellState<Scalar>& well_state)
const
1933 this->linSys_.extractCPRPressureMatrix(jacobian,
1944 template<
typename TypeTag>
1945 typename StandardWell<TypeTag>::EvalWell
1946 StandardWell<TypeTag>::
1947 pskinwater(
const Scalar throughput,
1948 const EvalWell& water_velocity,
1949 DeferredLogger& deferred_logger)
const
1951 if constexpr (Base::has_polymermw) {
1952 const int water_table_id = this->polymerWaterTable_();
1953 if (water_table_id <= 0) {
1954 OPM_DEFLOG_THROW(std::runtime_error,
1955 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1958 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1959 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1961 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1962 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1965 OPM_DEFLOG_THROW(std::runtime_error,
1966 fmt::format(
"Polymermw is not activated, while injecting "
1967 "skin pressure is requested for well {}", name()),
1976 template<
typename TypeTag>
1977 typename StandardWell<TypeTag>::EvalWell
1978 StandardWell<TypeTag>::
1979 pskin(
const Scalar throughput,
1980 const EvalWell& water_velocity,
1981 const EvalWell& poly_inj_conc,
1982 DeferredLogger& deferred_logger)
const
1984 if constexpr (Base::has_polymermw) {
1985 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1986 const EvalWell water_velocity_abs = abs(water_velocity);
1987 if (poly_inj_conc == 0.) {
1988 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1990 const int polymer_table_id = this->polymerTable_();
1991 if (polymer_table_id <= 0) {
1992 OPM_DEFLOG_THROW(std::runtime_error,
1993 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1996 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1997 const Scalar reference_concentration = skprpolytable.refConcentration;
1998 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2000 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2001 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2002 if (poly_inj_conc == reference_concentration) {
2003 return sign * pskin_poly;
2006 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2007 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2008 return sign * pskin;
2010 OPM_DEFLOG_THROW(std::runtime_error,
2011 fmt::format(
"Polymermw is not activated, while injecting "
2012 "skin pressure is requested for well {}", name()),
2021 template<
typename TypeTag>
2022 typename StandardWell<TypeTag>::EvalWell
2023 StandardWell<TypeTag>::
2024 wpolymermw(
const Scalar throughput,
2025 const EvalWell& water_velocity,
2026 DeferredLogger& deferred_logger)
const
2028 if constexpr (Base::has_polymermw) {
2029 const int table_id = this->polymerInjTable_();
2030 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2031 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2032 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2033 if (this->wpolymer() == 0.) {
2034 return molecular_weight;
2036 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2037 return molecular_weight;
2039 OPM_DEFLOG_THROW(std::runtime_error,
2040 fmt::format(
"Polymermw is not activated, while injecting "
2041 "polymer molecular weight is requested for well {}", name()),
2050 template<
typename TypeTag>
2052 StandardWell<TypeTag>::
2053 updateWaterThroughput([[maybe_unused]]
const double dt,
2054 WellState<Scalar>& well_state)
const
2056 if constexpr (Base::has_polymermw) {
2057 if (!this->isInjector()) {
2061 auto& perf_water_throughput = well_state.well(this->index_of_well_)
2062 .perf_data.water_throughput;
2064 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2065 const Scalar perf_water_vel =
2066 this->primary_variables_.value(Bhp + 1 + perf);
2070 if (perf_water_vel > Scalar{0}) {
2071 perf_water_throughput[perf] += perf_water_vel * dt;
2081 template<
typename TypeTag>
2083 StandardWell<TypeTag>::
2084 handleInjectivityRate(
const Simulator& simulator,
2086 std::vector<EvalWell>& cq_s)
const
2088 const int cell_idx = this->well_cells_[perf];
2089 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2090 const auto& fs = int_quants.fluidState();
2091 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2092 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2093 const int wat_vel_index = Bhp + 1 + perf;
2094 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2098 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2104 template<
typename TypeTag>
2106 StandardWell<TypeTag>::
2107 handleInjectivityEquations(
const Simulator& simulator,
2108 const WellState<Scalar>& well_state,
2110 const EvalWell& water_flux_s,
2111 DeferredLogger& deferred_logger)
2113 const int cell_idx = this->well_cells_[perf];
2114 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2115 const auto& fs = int_quants.fluidState();
2116 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2117 const EvalWell water_flux_r = water_flux_s / b_w;
2118 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2119 const EvalWell water_velocity = water_flux_r / area;
2120 const int wat_vel_index = Bhp + 1 + perf;
2123 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2125 const auto& ws = well_state.well(this->index_of_well_);
2126 const auto& perf_data = ws.perf_data;
2127 const auto& perf_water_throughput = perf_data.water_throughput;
2128 const Scalar throughput = perf_water_throughput[perf];
2129 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2131 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2132 poly_conc.setValue(this->wpolymer());
2135 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2136 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2138 StandardWellAssemble<FluidSystem,Indices>(*this).
2139 assembleInjectivityEq(eq_pskin,
2144 this->primary_variables_.numWellEq(),
2152 template<
typename TypeTag>
2154 StandardWell<TypeTag>::
2155 checkConvergenceExtraEqs(
const std::vector<Scalar>& res,
2156 ConvergenceReport& report)
const
2161 if constexpr (Base::has_polymermw) {
2162 WellConvergence(*this).
2163 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2171 template<
typename TypeTag>
2173 StandardWell<TypeTag>::
2174 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2175 const IntensiveQuantities& int_quants,
2176 const WellState<Scalar>& well_state,
2178 std::vector<RateVector>& connectionRates,
2179 DeferredLogger& deferred_logger)
const
2182 EvalWell cq_s_polymw = cq_s_poly;
2183 if (this->isInjector()) {
2184 const int wat_vel_index = Bhp + 1 + perf;
2185 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2186 if (water_velocity > 0.) {
2187 const auto& ws = well_state.well(this->index_of_well_);
2188 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2189 const Scalar throughput = perf_water_throughput[perf];
2190 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2191 cq_s_polymw *= molecular_weight;
2197 }
else if (this->isProducer()) {
2198 if (cq_s_polymw < 0.) {
2199 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2206 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2214 template<
typename TypeTag>
2215 std::optional<typename StandardWell<TypeTag>::Scalar>
2216 StandardWell<TypeTag>::
2217 computeBhpAtThpLimitProd(
const WellState<Scalar>& well_state,
2218 const Simulator& simulator,
2219 const SummaryState& summary_state,
2220 DeferredLogger& deferred_logger)
const
2222 return computeBhpAtThpLimitProdWithAlq(simulator,
2224 this->getALQ(well_state),
2229 template<
typename TypeTag>
2230 std::optional<typename StandardWell<TypeTag>::Scalar>
2231 StandardWell<TypeTag>::
2232 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2233 const SummaryState& summary_state,
2234 const Scalar alq_value,
2235 DeferredLogger& deferred_logger,
2236 bool iterate_if_no_solution)
const
2240 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2246 std::vector<Scalar> rates(3);
2247 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2248 this->adaptRatesForVFP(rates);
2252 Scalar max_pressure = 0.0;
2253 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2254 const int cell_idx = this->well_cells_[perf];
2255 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2256 const auto& fs = int_quants.fluidState();
2257 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2258 max_pressure = std::max(max_pressure, pressure_cell);
2260 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2263 this->connections_.rho(),
2265 this->getTHPConstraint(summary_state),
2269 auto v = frates(*bhpAtLimit);
2270 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2275 if (!iterate_if_no_solution)
2276 return std::nullopt;
2278 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2282 std::vector<Scalar> rates(3);
2283 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2284 this->adaptRatesForVFP(rates);
2288 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2291 this->connections_.rho(),
2293 this->getTHPConstraint(summary_state),
2299 auto v = frates(*bhpAtLimit);
2300 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2306 return std::nullopt;
2311 template<
typename TypeTag>
2312 std::optional<typename StandardWell<TypeTag>::Scalar>
2313 StandardWell<TypeTag>::
2314 computeBhpAtThpLimitInj(
const Simulator& simulator,
2315 const SummaryState& summary_state,
2316 DeferredLogger& deferred_logger)
const
2319 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2325 std::vector<Scalar> rates(3);
2326 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2330 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2332 this->connections_.rho(),
2343 template<
typename TypeTag>
2345 StandardWell<TypeTag>::
2346 iterateWellEqWithControl(
const Simulator& simulator,
2348 const Well::InjectionControls& inj_controls,
2349 const Well::ProductionControls& prod_controls,
2350 WellState<Scalar>& well_state,
2351 const GroupState<Scalar>& group_state,
2352 DeferredLogger& deferred_logger)
2354 const int max_iter = this->param_.max_inner_iter_wells_;
2357 bool relax_convergence =
false;
2358 this->regularize_ =
false;
2360 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2362 if (it > this->param_.strict_inner_iter_wells_) {
2363 relax_convergence =
true;
2364 this->regularize_ =
true;
2367 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2369 converged = report.converged();
2375 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2382 }
while (it < max_iter);
2388 template<
typename TypeTag>
2390 StandardWell<TypeTag>::
2391 iterateWellEqWithSwitching(
const Simulator& simulator,
2393 const Well::InjectionControls& inj_controls,
2394 const Well::ProductionControls& prod_controls,
2395 WellState<Scalar>& well_state,
2396 const GroupState<Scalar>& group_state,
2397 DeferredLogger& deferred_logger,
2398 const bool fixed_control ,
2399 const bool fixed_status )
2401 const int max_iter = this->param_.max_inner_iter_wells_;
2403 bool converged =
false;
2404 bool relax_convergence =
false;
2405 this->regularize_ =
false;
2406 const auto& summary_state = simulator.vanguard().summaryState();
2411 constexpr int min_its_after_switch = 4;
2412 int its_since_last_switch = min_its_after_switch;
2413 int switch_count= 0;
2415 const auto well_status_orig = this->wellStatus_;
2416 const auto operability_orig = this->operability_status_;
2417 auto well_status_cur = well_status_orig;
2418 int status_switch_count = 0;
2420 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2421 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2423 const bool allow_switching =
2424 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2425 (!fixed_control || !fixed_status) && allow_open;
2427 bool changed =
false;
2428 bool final_check =
false;
2430 this->operability_status_.resetOperability();
2431 this->operability_status_.solvable =
true;
2433 its_since_last_switch++;
2434 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2435 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2436 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2437 inj_controls, prod_controls, wqTotal,
2438 deferred_logger, fixed_control, fixed_status);
2440 its_since_last_switch = 0;
2442 if (well_status_cur != this->wellStatus_) {
2443 well_status_cur = this->wellStatus_;
2444 status_switch_count++;
2447 if (!changed && final_check) {
2450 final_check =
false;
2454 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2456 if (it > this->param_.strict_inner_iter_wells_) {
2457 relax_convergence =
true;
2458 this->regularize_ =
true;
2461 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2463 converged = report.converged();
2467 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2469 its_since_last_switch = min_its_after_switch;
2476 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2478 }
while (it < max_iter);
2481 if (allow_switching){
2483 const bool is_stopped = this->wellIsStopped();
2484 if (this->wellHasTHPConstraints(summary_state)){
2485 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2486 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2488 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2492 this->wellStatus_ = well_status_orig;
2493 this->operability_status_ = operability_orig;
2494 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2495 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2496 deferred_logger.debug(message);
2502 template<
typename TypeTag>
2503 std::vector<typename StandardWell<TypeTag>::Scalar>
2509 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2510 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2511 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2512 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2513 const int cell_idx = this->well_cells_[perf];
2514 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2515 std::vector<Scalar> mob(this->num_components_, 0.);
2516 getMobility(simulator, perf, mob, deferred_logger);
2517 std::vector<Scalar> cq_s(this->num_components_, 0.);
2518 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2519 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2520 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2522 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2523 cq_s, perf_rates, deferred_logger);
2524 for (
int comp = 0; comp < this->num_components_; ++comp) {
2525 well_q_s[comp] += cq_s[comp];
2528 const auto& comm = this->parallel_well_info_.communication();
2529 if (comm.size() > 1)
2531 comm.sum(well_q_s.data(), well_q_s.size());
2538 template <
typename TypeTag>
2539 std::vector<typename StandardWell<TypeTag>::Scalar>
2543 const int num_pri_vars = this->primary_variables_.numWellEq();
2544 std::vector<Scalar> retval(num_pri_vars);
2545 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2546 retval[ii] = this->primary_variables_.value(ii);
2555 template <
typename TypeTag>
2557 StandardWell<TypeTag>::
2558 setPrimaryVars(
typename std::vector<Scalar>::const_iterator it)
2560 const int num_pri_vars = this->primary_variables_.numWellEq();
2561 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2562 this->primary_variables_.setValue(ii, it[ii]);
2564 return num_pri_vars;
2568 template <
typename TypeTag>
2569 typename StandardWell<TypeTag>::Eval
2570 StandardWell<TypeTag>::
2571 connectionRateEnergy(
const Scalar maxOilSaturation,
2572 const std::vector<EvalWell>& cq_s,
2573 const IntensiveQuantities& intQuants,
2574 DeferredLogger& deferred_logger)
const
2576 auto fs = intQuants.fluidState();
2578 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2579 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2584 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2585 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2586 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2587 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2588 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2591 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2592 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2597 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2599 deferred_logger.debug(
2600 fmt::format(
"Problematic d value {} obtained for well {}"
2601 " during calculateSinglePerf with rs {}"
2602 ", rv {}. Continue as if no dissolution (rs = 0) and"
2603 " vaporization (rv = 0) for this connection.",
2604 d, this->name(), fs.Rs(), fs.Rv()));
2605 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2607 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2608 cq_r_thermal = (cq_s[gasCompIdx] -
2609 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2610 (d * this->extendEval(fs.invB(phaseIdx)) );
2611 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2613 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2615 (d * this->extendEval(fs.invB(phaseIdx)) );
2621 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2623 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2624 fs.setTemperature(this->well_ecl_.inj_temperature());
2625 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2626 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2627 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2628 paramCache.setRegionIndex(pvtRegionIdx);
2629 paramCache.setMaxOilSat(maxOilSaturation);
2630 paramCache.updatePhase(fs, phaseIdx);
2632 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2633 fs.setDensity(phaseIdx, rho);
2634 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2635 fs.setEnthalpy(phaseIdx, h);
2636 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2637 result += getValue(cq_r_thermal);
2638 }
else if (cq_r_thermal > 0.0) {
2639 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2640 result += Base::restrictEval(cq_r_thermal);
2643 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2644 result += Base::restrictEval(cq_r_thermal);
2648 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:97
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:312
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:299
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:1424
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