23#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
24#define OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
26#include <opm/simulators/wells/StandardWell.hpp>
29#include <opm/common/Exceptions.hpp>
31#include <opm/input/eclipse/Units/Units.hpp>
33#include <opm/material/densead/EvaluationFormat.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>
41#include <fmt/format.h>
50template<
class dValue,
class Value>
51auto dValueError(
const dValue& d,
52 const std::string& name,
53 const std::string& methodName,
56 const Value& pressure)
58 return fmt::format(
"Problematic d value {} obtained for well {}"
59 " during {} calculations with rs {}"
60 ", rv {} and pressure {}."
61 " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
62 " for this connection.", d, name, methodName, Rs, Rv, pressure);
70 template<
typename TypeTag>
71 StandardWell<TypeTag>::
72 StandardWell(
const Well& well,
73 const ParallelWellInfo<Scalar>& pw_info,
75 const ModelParameters& param,
76 const RateConverterType& rate_converter,
77 const int pvtRegionIdx,
78 const int num_components,
80 const int index_of_well,
81 const std::vector<PerforationData<Scalar>>& perf_data)
82 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
83 , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices>&>(*this))
86 assert(this->num_components_ == numWellConservationEq);
93 template<
typename TypeTag>
95 StandardWell<TypeTag>::
96 init(
const PhaseUsage* phase_usage_arg,
97 const std::vector<Scalar>& depth_arg,
98 const Scalar gravity_arg,
99 const std::vector< Scalar >& B_avg,
100 const bool changed_to_open_this_step)
102 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
103 this->StdWellEval::init(this->perf_depth_, depth_arg, Base::has_polymermw);
110 template<
typename TypeTag>
111 void StandardWell<TypeTag>::
112 initPrimaryVariablesEvaluation()
114 this->primary_variables_.init();
121 template<
typename TypeTag>
122 template<
class Value>
124 StandardWell<TypeTag>::
125 computePerfRate(
const IntensiveQuantities& intQuants,
126 const std::vector<Value>& mob,
128 const std::vector<Scalar>& Tw,
131 std::vector<Value>& cq_s,
132 PerforationRates<Scalar>& perf_rates,
133 DeferredLogger& deferred_logger)
const
135 auto obtain = [
this](
const Eval& value)
137 if constexpr (std::is_same_v<Value, Scalar>) {
138 static_cast<void>(
this);
139 return getValue(value);
141 return this->extendEval(value);
144 auto obtainN = [](
const auto& value)
146 if constexpr (std::is_same_v<Value, Scalar>) {
147 return getValue(value);
152 auto zeroElem = [
this]()
154 if constexpr (std::is_same_v<Value, Scalar>) {
155 static_cast<void>(
this);
158 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
162 const auto& fs = intQuants.fluidState();
163 const Value pressure = obtain(this->getPerfCellPressure(fs));
164 const Value rs = obtain(fs.Rs());
165 const Value rv = obtain(fs.Rv());
166 const Value rvw = obtain(fs.Rvw());
167 const Value rsw = obtain(fs.Rsw());
169 std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
170 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
171 if (!FluidSystem::phaseIsActive(phaseIdx)) {
174 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
175 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
177 if constexpr (has_solvent) {
178 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
181 if constexpr (has_zFraction) {
182 if (this->isInjector()) {
183 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
184 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
185 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
189 Value skin_pressure = zeroElem();
191 if (this->isInjector()) {
192 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
193 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
198 std::vector<Value> cmix_s(this->numComponents(), zeroElem());
199 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
200 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
221 template<
typename TypeTag>
222 template<
class Value>
224 StandardWell<TypeTag>::
225 computePerfRate(
const std::vector<Value>& mob,
226 const Value& pressure,
232 std::vector<Value>& b_perfcells_dense,
233 const std::vector<Scalar>& Tw,
236 const Value& skin_pressure,
237 const std::vector<Value>& cmix_s,
238 std::vector<Value>& cq_s,
239 PerforationRates<Scalar>& perf_rates,
240 DeferredLogger& deferred_logger)
const
243 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
244 Value drawdown = pressure - well_pressure;
245 if (this->isInjector()) {
246 drawdown += skin_pressure;
252 if (!allow_cf && this->isInjector()) {
257 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
258 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
259 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
262 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
263 gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw);
264 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
265 gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw);
269 if (!allow_cf && this->isProducer()) {
274 Value total_mob_dense = mob[0];
275 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
276 total_mob_dense += mob[componentIdx];
280 Value volumeRatio = bhp * 0.0;
282 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
283 disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
284 cmix_s, b_perfcells_dense, deferred_logger);
288 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
289 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
290 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
293 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
294 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
295 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
298 if constexpr (Indices::enableSolvent) {
299 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
302 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
303 gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
304 cmix_s, b_perfcells_dense, deferred_logger);
306 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
307 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
308 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
310 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
311 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
312 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
318 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
319 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
320 Value cqt_is = cqt_i / volumeRatio;
321 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
325 if (this->isProducer()) {
326 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
327 gasOilPerfRateInj(cq_s, perf_rates,
328 rv, rs, pressure, rvw, deferred_logger);
330 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
332 gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
333 pressure, deferred_logger);
340 template<
typename TypeTag>
342 StandardWell<TypeTag>::
343 assembleWellEqWithoutIteration(
const Simulator& simulator,
345 const Well::InjectionControls& inj_controls,
346 const Well::ProductionControls& prod_controls,
347 WellState<Scalar>& well_state,
348 const GroupState<Scalar>& group_state,
349 DeferredLogger& deferred_logger)
353 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
356 this->linSys_.clear();
358 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
359 prod_controls, well_state,
360 group_state, deferred_logger);
366 template<
typename TypeTag>
368 StandardWell<TypeTag>::
369 assembleWellEqWithoutIterationImpl(
const Simulator& simulator,
371 const Well::InjectionControls& inj_controls,
372 const Well::ProductionControls& prod_controls,
373 WellState<Scalar>& well_state,
374 const GroupState<Scalar>& group_state,
375 DeferredLogger& deferred_logger)
378 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
379 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
381 auto& ws = well_state.well(this->index_of_well_);
382 ws.phase_mixing_rates.fill(0.0);
385 const int np = this->number_of_phases_;
387 std::vector<RateVector> connectionRates = this->connectionRates_;
389 auto& perf_data = ws.perf_data;
390 auto& perf_rates = perf_data.phase_rates;
391 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
393 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
394 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
395 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
396 calculateSinglePerf(simulator, perf, well_state, connectionRates,
397 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
400 if constexpr (has_polymer && Base::has_polymermw) {
401 if (this->isInjector()) {
402 handleInjectivityEquations(simulator, well_state, perf,
403 water_flux_s, deferred_logger);
406 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
408 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
410 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
412 StandardWellAssemble<FluidSystem,Indices>(*this).
413 assemblePerforationEq(cq_s_effective,
416 this->primary_variables_.numWellEq(),
420 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
421 auto& perf_rate_solvent = perf_data.solvent_rates;
422 perf_rate_solvent[perf] = cq_s[componentIdx].value();
424 perf_rates[perf*np + this->modelCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
428 if constexpr (has_zFraction) {
429 StandardWellAssemble<FluidSystem,Indices>(*this).
430 assembleZFracEq(cq_s_zfrac_effective,
432 this->primary_variables_.numWellEq(),
437 this->connectionRates_ = connectionRates;
442 const auto& comm = this->parallel_well_info_.communication();
443 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
447 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
450 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
453 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
454 if (FluidSystem::numActivePhases() > 1) {
456 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
457 this->F0_[componentIdx]) * volume / dt;
459 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
460 StandardWellAssemble<FluidSystem,Indices>(*this).
461 assembleSourceEq(resWell_loc,
463 this->primary_variables_.numWellEq(),
467 const auto& summaryState = simulator.vanguard().summaryState();
468 const Schedule& schedule = simulator.vanguard().schedule();
469 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
470 StandardWellAssemble<FluidSystem,Indices>(*this).
471 assembleControlEq(well_state, group_state,
472 schedule, summaryState,
473 inj_controls, prod_controls,
474 this->primary_variables_,
475 this->connections_.rho(),
477 stopped_or_zero_target,
483 this->linSys_.invert();
485 OPM_DEFLOG_PROBLEM(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
492 template<
typename TypeTag>
494 StandardWell<TypeTag>::
495 calculateSinglePerf(
const Simulator& simulator,
497 WellState<Scalar>& well_state,
498 std::vector<RateVector>& connectionRates,
499 std::vector<EvalWell>& cq_s,
500 EvalWell& water_flux_s,
501 EvalWell& cq_s_zfrac_effective,
502 DeferredLogger& deferred_logger)
const
504 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
505 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
506 const int cell_idx = this->well_cells_[perf];
507 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
508 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
509 getMobility(simulator, perf, mob, deferred_logger);
511 PerforationRates<Scalar> perf_rates;
512 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
513 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
514 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
515 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
516 cq_s, perf_rates, deferred_logger);
518 auto& ws = well_state.well(this->index_of_well_);
519 auto& perf_data = ws.perf_data;
520 if constexpr (has_polymer && Base::has_polymermw) {
521 if (this->isInjector()) {
524 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
525 water_flux_s = cq_s[water_comp_idx];
528 handleInjectivityRate(simulator, perf, cq_s);
533 if (this->isProducer()) {
534 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
535 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
536 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
537 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
538 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.dis_gas;
539 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.dis_gas_in_water;
540 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.vap_oil;
541 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.vap_wat;
544 if constexpr (has_energy) {
545 connectionRates[perf][Indices::contiEnergyEqIdx] =
546 connectionRateEnergy(simulator.problem().maxOilSaturation(cell_idx),
547 cq_s, intQuants, deferred_logger);
550 if constexpr (has_polymer) {
551 std::variant<Scalar,EvalWell> polymerConcentration;
552 if (this->isInjector()) {
553 polymerConcentration = this->wpolymer();
555 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
556 intQuants.polymerViscosityCorrection());
559 [[maybe_unused]] EvalWell cq_s_poly;
560 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
562 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
563 cq_s, polymerConcentration);
565 if constexpr (Base::has_polymermw) {
566 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
567 perf, connectionRates, deferred_logger);
571 if constexpr (has_foam) {
572 std::variant<Scalar,EvalWell> foamConcentration;
573 if (this->isInjector()) {
574 foamConcentration = this->wfoam();
576 foamConcentration = this->extendEval(intQuants.foamConcentration());
578 connectionRates[perf][Indices::contiFoamEqIdx] =
579 this->connections_.connectionRateFoam(cq_s, foamConcentration,
580 FoamModule::transportPhase(),
584 if constexpr (has_zFraction) {
585 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
586 if (this->isInjector()) {
587 solventConcentration = this->wsolvent();
589 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
590 this->extendEval(intQuants.yVolume())};
592 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
593 cq_s_zfrac_effective) =
594 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
595 perf_rates.dis_gas, cq_s,
596 solventConcentration);
599 if constexpr (has_brine) {
600 std::variant<Scalar,EvalWell> saltConcentration;
601 if (this->isInjector()) {
602 saltConcentration = this->wsalt();
604 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
607 connectionRates[perf][Indices::contiBrineEqIdx] =
608 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
609 perf_rates.vap_wat, cq_s,
613 if constexpr (has_micp) {
614 std::variant<Scalar,EvalWell> microbialConcentration;
615 std::variant<Scalar,EvalWell> oxygenConcentration;
616 std::variant<Scalar,EvalWell> ureaConcentration;
617 if (this->isInjector()) {
618 microbialConcentration = this->wmicrobes();
619 oxygenConcentration = this->woxygen();
620 ureaConcentration = this->wurea();
622 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
623 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
624 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
626 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
627 connectionRates[perf][Indices::contiOxygenEqIdx],
628 connectionRates[perf][Indices::contiUreaEqIdx]) =
629 this->connections_.connectionRatesMICP(cq_s,
630 microbialConcentration,
636 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
641 template<
typename TypeTag>
642 template<
class Value>
644 StandardWell<TypeTag>::
645 getMobility(
const Simulator& simulator,
647 std::vector<Value>& mob,
648 DeferredLogger& deferred_logger)
const
650 auto obtain = [
this](
const Eval& value)
652 if constexpr (std::is_same_v<Value, Scalar>) {
653 static_cast<void>(
this);
654 return getValue(value);
656 return this->extendEval(value);
659 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
660 obtain, deferred_logger);
663 if constexpr (has_polymer) {
664 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
665 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
670 if constexpr (!Base::has_polymermw) {
671 if constexpr (std::is_same_v<Value, Scalar>) {
672 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
673 for (std::size_t i = 0; i < mob.size(); ++i) {
674 mob_eval[i].setValue(mob[i]);
676 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
677 for (std::size_t i = 0; i < mob.size(); ++i) {
678 mob[i] = getValue(mob_eval[i]);
681 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
687 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
688 const Scalar bhp = this->primary_variables_.value(Bhp);
689 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
690 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
691 for (std::size_t i = 0; i < mob.size(); ++i) {
692 mob[i] *= multiplier;
698 template<
typename TypeTag>
700 StandardWell<TypeTag>::
701 updateWellState(
const Simulator& simulator,
702 const BVectorWell& dwells,
703 WellState<Scalar>& well_state,
704 DeferredLogger& deferred_logger)
706 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
708 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
709 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
711 const auto& summary_state = simulator.vanguard().summaryState();
712 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
713 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
720 template<
typename TypeTag>
722 StandardWell<TypeTag>::
723 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
724 const bool stop_or_zero_rate_target,
725 DeferredLogger& deferred_logger)
727 const Scalar dFLimit = this->param_.dwell_fraction_max_;
728 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
729 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
732 if constexpr (Base::has_polymermw) {
733 this->primary_variables_.updateNewtonPolyMW(dwells);
736 this->primary_variables_.checkFinite(deferred_logger);
743 template<
typename TypeTag>
745 StandardWell<TypeTag>::
746 updateWellStateFromPrimaryVariables(WellState<Scalar>& well_state,
747 const SummaryState& summary_state,
748 DeferredLogger& deferred_logger)
const
750 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
753 if constexpr (Base::has_polymermw) {
754 this->primary_variables_.copyToWellStatePolyMW(well_state);
762 template<
typename TypeTag>
764 StandardWell<TypeTag>::
765 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
770 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
771 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
773 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
774 std::vector<Scalar> mob(this->num_components_, 0.0);
775 getMobility(simulator, perf, mob, deferred_logger);
777 const int cell_idx = this->well_cells_[perf];
778 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
779 const auto& fs = int_quantities.fluidState();
781 Scalar p_r = this->getPerfCellPressure(fs).value();
784 std::vector<Scalar> b_perf(this->num_components_);
785 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
786 if (!FluidSystem::phaseIsActive(phase)) {
789 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
790 b_perf[comp_idx] = fs.invB(phase).value();
792 if constexpr (has_solvent) {
793 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
797 const Scalar h_perf = this->connections_.pressure_diff(perf);
798 const Scalar pressure_diff = p_r - h_perf;
803 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
804 deferred_logger.debug(
"CROSSFLOW_IPR",
805 "cross flow found when updateIPR for well " + name()
806 +
" . The connection is ignored in IPR calculations");
812 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
813 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
814 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
818 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
819 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
820 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
821 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
822 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
823 ipr_b_perf[comp_idx] += tw_mob;
827 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
828 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
829 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
830 const Scalar rs = (fs.Rs()).value();
831 const Scalar rv = (fs.Rv()).value();
833 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
834 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
836 ipr_a_perf[gas_comp_idx] += dis_gas_a;
837 ipr_a_perf[oil_comp_idx] += vap_oil_a;
839 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
840 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
842 ipr_b_perf[gas_comp_idx] += dis_gas_b;
843 ipr_b_perf[oil_comp_idx] += vap_oil_b;
846 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
847 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
848 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
851 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
852 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
855 template<
typename TypeTag>
857 StandardWell<TypeTag>::
858 updateIPRImplicit(
const Simulator& simulator,
859 WellState<Scalar>& well_state,
860 DeferredLogger& deferred_logger)
869 auto rates = well_state.well(this->index_of_well_).surface_rates;
871 for (std::size_t p = 0; p < rates.size(); ++p) {
872 zero_rates &= rates[p] == 0.0;
874 auto& ws = well_state.well(this->index_of_well_);
876 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
877 deferred_logger.debug(msg);
889 const auto& group_state = simulator.problem().wellModel().groupState();
891 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
892 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
894 auto inj_controls = Well::InjectionControls(0);
895 auto prod_controls = Well::ProductionControls(0);
896 prod_controls.addControl(Well::ProducerCMode::BHP);
897 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
900 const auto cmode = ws.production_cmode;
901 ws.production_cmode = Well::ProducerCMode::BHP;
902 const double dt = simulator.timeStepSize();
903 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
905 const size_t nEq = this->primary_variables_.numWellEq();
909 for (
size_t i=0; i < nEq; ++i){
914 BVectorWell x_well(1);
915 x_well[0].resize(nEq);
916 this->linSys_.solve(rhs, x_well);
918 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
919 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
920 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
921 for (
size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
923 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
925 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
928 ws.production_cmode = cmode;
931 template<
typename TypeTag>
933 StandardWell<TypeTag>::
934 checkOperabilityUnderBHPLimit(
const WellState<Scalar>& well_state,
935 const Simulator& simulator,
936 DeferredLogger& deferred_logger)
938 const auto& summaryState = simulator.vanguard().summaryState();
939 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
942 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
943 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
946 Scalar total_ipr_mass_rate = 0.0;
947 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
949 if (!FluidSystem::phaseIsActive(phaseIdx)) {
953 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
954 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
956 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
957 total_ipr_mass_rate += ipr_rate * rho;
959 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
960 this->operability_status_.operable_under_only_bhp_limit =
false;
964 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
968 std::vector<Scalar> well_rates_bhp_limit;
969 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
971 this->adaptRatesForVFP(well_rates_bhp_limit);
972 const Scalar thp_limit = this->getTHPConstraint(summaryState);
973 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
975 this->connections_.rho(),
976 this->getALQ(well_state),
979 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
980 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
991 this->operability_status_.operable_under_only_bhp_limit =
true;
992 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1000 template<
typename TypeTag>
1002 StandardWell<TypeTag>::
1003 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
1004 const WellState<Scalar>& well_state,
1005 DeferredLogger& deferred_logger)
1007 const auto& summaryState = simulator.vanguard().summaryState();
1008 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1009 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1012 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1014 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1015 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1016 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1018 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1019 if (this->isProducer() && *obtain_bhp < thp_limit) {
1020 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1021 +
" bars is SMALLER than thp limit "
1022 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1023 +
" bars as a producer for well " + name();
1024 deferred_logger.debug(msg);
1026 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1027 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1028 +
" bars is LARGER than thp limit "
1029 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1030 +
" bars as a injector for well " + name();
1031 deferred_logger.debug(msg);
1034 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1035 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1036 if (!this->wellIsStopped()) {
1037 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1038 deferred_logger.debug(
" could not find bhp value at thp limit "
1039 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1040 +
" bar for well " + name() +
", the well might need to be closed ");
1049 template<
typename TypeTag>
1051 StandardWell<TypeTag>::
1052 allDrawDownWrongDirection(
const Simulator& simulator)
const
1054 bool all_drawdown_wrong_direction =
true;
1056 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1057 const int cell_idx = this->well_cells_[perf];
1058 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1059 const auto& fs = intQuants.fluidState();
1061 const Scalar pressure = this->getPerfCellPressure(fs).value();
1062 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1065 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1066 const Scalar drawdown = pressure - well_pressure;
1071 if ( (drawdown < 0. && this->isInjector()) ||
1072 (drawdown > 0. && this->isProducer()) ) {
1073 all_drawdown_wrong_direction =
false;
1078 const auto& comm = this->parallel_well_info_.communication();
1079 if (comm.size() > 1)
1081 all_drawdown_wrong_direction =
1082 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1085 return all_drawdown_wrong_direction;
1091 template<
typename TypeTag>
1093 StandardWell<TypeTag>::
1094 canProduceInjectWithCurrentBhp(
const Simulator& simulator,
1095 const WellState<Scalar>& well_state,
1096 DeferredLogger& deferred_logger)
1098 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
1099 std::vector<Scalar> well_rates;
1100 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1102 const Scalar sign = (this->isProducer()) ? -1. : 1.;
1103 const Scalar threshold = sign * std::numeric_limits<Scalar>::min();
1105 bool can_produce_inject =
false;
1106 for (
const auto value : well_rates) {
1107 if (this->isProducer() && value < threshold) {
1108 can_produce_inject =
true;
1110 }
else if (this->isInjector() && value > threshold) {
1111 can_produce_inject =
true;
1116 if (!can_produce_inject) {
1117 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1120 return can_produce_inject;
1127 template<
typename TypeTag>
1129 StandardWell<TypeTag>::
1130 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1132 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1138 template<
typename TypeTag>
1139 typename StandardWell<TypeTag>::WellConnectionProps
1140 StandardWell<TypeTag>::
1141 computePropertiesForWellConnectionPressures(
const Simulator& simulator,
1142 const WellState<Scalar>& well_state)
const
1144 auto prop_func =
typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1146 [&model = simulator.model()](
int cell_idx,
int phase_idx)
1148 return model.intensiveQuantities(cell_idx, 0)
1149 .fluidState().temperature(phase_idx).value();
1153 [&model = simulator.model()](
int cell_idx)
1155 return model.intensiveQuantities(cell_idx, 0)
1156 .fluidState().saltConcentration().value();
1160 [&model = simulator.model()](
int cell_idx)
1162 return model.intensiveQuantities(cell_idx, 0)
1163 .fluidState().pvtRegionIndex();
1167 if constexpr (Indices::enableSolvent) {
1168 prop_func.solventInverseFormationVolumeFactor =
1169 [&model = simulator.model()](
int cell_idx)
1171 return model.intensiveQuantities(cell_idx, 0)
1172 .solventInverseFormationVolumeFactor().value();
1175 prop_func.solventRefDensity = [&model = simulator.model()](
int cell_idx)
1177 return model.intensiveQuantities(cell_idx, 0)
1178 .solventRefDensity();
1182 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1189 template<
typename TypeTag>
1194 const std::vector<Scalar>& B_avg,
1196 const bool relax_tolerance)
const
1200 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1202 Scalar tol_wells = this->param_.tolerance_wells_;
1204 constexpr Scalar stopped_factor = 1.e-4;
1206 constexpr Scalar dynamic_thp_factor = 1.e-1;
1207 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1208 tol_wells = tol_wells*stopped_factor;
1209 }
else if (this->getDynamicThpLimit()) {
1210 tol_wells = tol_wells*dynamic_thp_factor;
1213 std::vector<Scalar> res;
1216 this->param_.max_residual_allowed_,
1218 this->param_.relaxed_tolerance_flow_well_,
1220 this->wellIsStopped(),
1224 checkConvergenceExtraEqs(res, report);
1233 template<
typename TypeTag>
1241 auto fluidState = [&simulator,
this](
const int perf)
1243 const auto cell_idx = this->well_cells_[perf];
1244 return simulator.
model()
1245 .intensiveQuantities(cell_idx, 0).fluidState();
1248 const int np = this->number_of_phases_;
1249 auto setToZero = [np](Scalar* x) ->
void
1251 std::fill_n(x, np, 0.0);
1254 auto addVector = [np](
const Scalar* src, Scalar* dest) ->
void
1256 std::transform(src, src + np, dest, dest, std::plus<>{});
1259 auto& ws = well_state.well(this->index_of_well_);
1260 auto& perf_data = ws.perf_data;
1261 auto* wellPI = ws.productivity_index.data();
1262 auto* connPI = perf_data.prod_index.data();
1266 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1267 auto subsetPerfID = 0;
1269 for (
const auto& perf : *this->perf_data_) {
1270 auto allPerfID = perf.ecl_index;
1272 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) -> Scalar
1277 std::vector<Scalar> mob(this->num_components_, 0.0);
1278 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1280 const auto& fs = fluidState(subsetPerfID);
1283 if (this->isInjector()) {
1284 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1285 mob, connPI, deferred_logger);
1288 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1291 addVector(connPI, wellPI);
1298 const auto& comm = this->parallel_well_info_.communication();
1299 if (comm.size() > 1) {
1300 comm.sum(wellPI, np);
1303 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1304 "Internal logic error in processing connections for PI/II");
1309 template<
typename TypeTag>
1310 void StandardWell<TypeTag>::
1311 computeWellConnectionDensitesPressures(
const Simulator& simulator,
1312 const WellState<Scalar>& well_state,
1313 const WellConnectionProps& props,
1314 DeferredLogger& deferred_logger)
1319 const auto prop_func =
typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1324 [&model = simulator.model()](
const int cell,
1325 const std::vector<int>& phases,
1326 std::vector<Scalar>& mob)
1328 const auto& iq = model.intensiveQuantities(cell, 0);
1330 std::transform(phases.begin(), phases.end(), mob.begin(),
1331 [&iq](
const int phase) { return iq.mobility(phase).value(); });
1336 [&model = simulator.model()](
const int cell,
1337 const std::vector<int>& phases,
1338 std::vector<Scalar>& rho)
1340 const auto& fs = model.intensiveQuantities(cell, 0).fluidState();
1342 std::transform(phases.begin(), phases.end(), rho.begin(),
1343 [&fs](
const int phase) { return fs.density(phase).value(); });
1347 const auto stopped_or_zero_rate_target = this->
1348 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1351 .computeProperties(stopped_or_zero_rate_target, well_state,
1352 prop_func, props, deferred_logger);
1359 template<
typename TypeTag>
1361 StandardWell<TypeTag>::
1362 computeWellConnectionPressures(
const Simulator& simulator,
1363 const WellState<Scalar>& well_state,
1364 DeferredLogger& deferred_logger)
1366 const auto props = computePropertiesForWellConnectionPressures
1367 (simulator, well_state);
1369 computeWellConnectionDensitesPressures(simulator, well_state,
1370 props, deferred_logger);
1377 template<
typename TypeTag>
1379 StandardWell<TypeTag>::
1380 solveEqAndUpdateWellState(
const Simulator& simulator,
1381 WellState<Scalar>& well_state,
1382 DeferredLogger& deferred_logger)
1384 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1388 BVectorWell dx_well(1);
1389 dx_well[0].resize(this->primary_variables_.numWellEq());
1390 this->linSys_.solve( dx_well);
1392 updateWellState(simulator, dx_well, well_state, deferred_logger);
1399 template<
typename TypeTag>
1401 StandardWell<TypeTag>::
1402 calculateExplicitQuantities(
const Simulator& simulator,
1403 const WellState<Scalar>& well_state,
1404 DeferredLogger& deferred_logger)
1406 updatePrimaryVariables(simulator, well_state, deferred_logger);
1407 initPrimaryVariablesEvaluation();
1408 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1409 this->computeAccumWell();
1414 template<
typename TypeTag>
1417 apply(
const BVector& x, BVector& Ax)
const
1419 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1421 if (this->param_.matrix_add_well_contributions_)
1427 this->linSys_.
apply(x, Ax);
1433 template<
typename TypeTag>
1436 apply(BVector& r)
const
1438 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1440 this->linSys_.
apply(r);
1446 template<
typename TypeTag>
1454 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1457 xw[0].resize(this->primary_variables_.numWellEq());
1459 this->linSys_.recoverSolutionWell(x, xw);
1460 updateWellState(simulator, xw, well_state, deferred_logger);
1466 template<
typename TypeTag>
1471 std::vector<Scalar>& well_flux,
1475 const int np = this->number_of_phases_;
1476 well_flux.resize(np, 0.0);
1478 const bool allow_cf = this->getAllowCrossFlow();
1480 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1481 const int cell_idx = this->well_cells_[perf];
1482 const auto& intQuants = simulator.
model().intensiveQuantities(cell_idx, 0);
1484 std::vector<Scalar> mob(this->num_components_, 0.);
1485 getMobility(simulator, perf, mob, deferred_logger);
1486 Scalar trans_mult = simulator.
problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1487 const auto& wellstate_nupcol = simulator.
problem().wellModel().nupcolWellState().well(this->index_of_well_);
1488 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1490 std::vector<Scalar> cq_s(this->num_components_, 0.);
1492 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1493 cq_s, perf_rates, deferred_logger);
1495 for(
int p = 0; p < np; ++p) {
1496 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1500 if constexpr (has_solvent) {
1501 const auto& pu = this->phaseUsage();
1502 assert(pu.phase_used[Gas]);
1503 const int gas_pos = pu.phase_pos[Gas];
1504 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1507 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1512 template<
typename TypeTag>
1514 StandardWell<TypeTag>::
1515 computeWellRatesWithBhpIterations(
const Simulator& simulator,
1517 std::vector<Scalar>& well_flux,
1518 DeferredLogger& deferred_logger)
const
1522 StandardWell<TypeTag> well_copy(*
this);
1523 well_copy.resetDampening();
1528 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1529 const auto& group_state = simulator.problem().wellModel().groupState();
1532 const auto& summary_state = simulator.vanguard().summaryState();
1533 auto inj_controls = well_copy.well_ecl_.isInjector()
1534 ? well_copy.well_ecl_.injectionControls(summary_state)
1535 : Well::InjectionControls(0);
1536 auto prod_controls = well_copy.well_ecl_.isProducer()
1537 ? well_copy.well_ecl_.productionControls(summary_state) :
1538 Well::ProductionControls(0);
1541 auto& ws = well_state_copy.well(this->index_of_well_);
1542 if (well_copy.well_ecl_.isInjector()) {
1543 inj_controls.bhp_limit = bhp;
1544 ws.injection_cmode = Well::InjectorCMode::BHP;
1546 prod_controls.bhp_limit = bhp;
1547 ws.production_cmode = Well::ProducerCMode::BHP;
1552 const int np = this->number_of_phases_;
1553 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1554 for (
int phase = 0; phase < np; ++phase){
1555 well_state_copy.wellRates(this->index_of_well_)[phase]
1556 = sign * ws.well_potentials[phase];
1558 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1559 well_copy.initPrimaryVariablesEvaluation();
1560 well_copy.computeAccumWell();
1562 const double dt = simulator.timeStepSize();
1563 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1565 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1566 " potentials are computed based on unconverged solution";
1567 deferred_logger.debug(msg);
1569 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1570 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1571 well_copy.initPrimaryVariablesEvaluation();
1572 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1578 template<
typename TypeTag>
1579 std::vector<typename StandardWell<TypeTag>::Scalar>
1580 StandardWell<TypeTag>::
1581 computeWellPotentialWithTHP(
const Simulator& simulator,
1582 DeferredLogger& deferred_logger,
1583 const WellState<Scalar>& well_state)
const
1585 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1586 const auto& summary_state = simulator.vanguard().summaryState();
1588 const auto& well = this->well_ecl_;
1589 if (well.isInjector()){
1590 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1591 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1592 if (bhp_at_thp_limit) {
1593 const Scalar bhp = std::min(*bhp_at_thp_limit,
1594 static_cast<Scalar
>(controls.bhp_limit));
1595 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1597 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1598 "Failed in getting converged thp based potential calculation for well "
1599 + name() +
". Instead the bhp based value is used");
1600 const Scalar bhp = controls.bhp_limit;
1601 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1604 computeWellRatesWithThpAlqProd(
1605 simulator, summary_state,
1606 deferred_logger, potentials, this->getALQ(well_state)
1613 template<
typename TypeTag>
1615 StandardWell<TypeTag>::
1616 computeWellPotentialsImplicit(
const Simulator& simulator,
1617 std::vector<Scalar>& well_potentials,
1618 DeferredLogger& deferred_logger)
const
1623 StandardWell<TypeTag> well_copy(*
this);
1626 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1627 const auto& group_state = simulator.problem().wellModel().groupState();
1628 auto& ws = well_state_copy.well(this->index_of_well_);
1631 const auto& summary_state = simulator.vanguard().summaryState();
1632 auto inj_controls = well_copy.well_ecl_.isInjector()
1633 ? well_copy.well_ecl_.injectionControls(summary_state)
1634 : Well::InjectionControls(0);
1635 auto prod_controls = well_copy.well_ecl_.isProducer()
1636 ? well_copy.well_ecl_.productionControls(summary_state) :
1637 Well::ProductionControls(0);
1640 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1643 const int num_perf = ws.perf_data.size();
1644 for (
int perf = 0; perf < num_perf; ++perf) {
1645 ws.perf_data.pressure[perf] = ws.bhp + well_copy.connections_.pressure_diff(perf);
1648 const int np = this->number_of_phases_;
1649 bool trivial =
true;
1650 for (
int phase = 0; phase < np; ++phase){
1651 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1654 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1655 for (
int phase = 0; phase < np; ++phase) {
1656 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1660 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1661 const double dt = simulator.timeStepSize();
1663 bool converged =
false;
1664 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1665 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1667 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1671 well_potentials.clear();
1672 well_potentials.resize(np, 0.0);
1673 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1674 if (has_solvent && comp_idx == Indices::contiSolventEqIdx)
continue;
1675 const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
1676 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1680 if constexpr (has_solvent) {
1682 assert(pu.phase_used[Gas]);
1683 const int gas_pos = pu.phase_pos[Gas];
1684 const EvalWell rate = well_copy.primary_variables_.getQs(Indices::contiSolventEqIdx);
1685 well_potentials[gas_pos] += rate.value();
1691 template<
typename TypeTag>
1692 typename StandardWell<TypeTag>::Scalar
1693 StandardWell<TypeTag>::
1694 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &simulator,
1695 const SummaryState &summary_state,
1696 DeferredLogger& deferred_logger,
1697 std::vector<Scalar>& potentials,
1701 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1702 simulator, summary_state, alq, deferred_logger);
1703 if (bhp_at_thp_limit) {
1704 const auto& controls = this->well_ecl_.productionControls(summary_state);
1705 bhp = std::max(*bhp_at_thp_limit,
1706 static_cast<Scalar
>(controls.bhp_limit));
1707 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1710 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1711 "Failed in getting converged thp based potential calculation for well "
1712 + name() +
". Instead the bhp based value is used");
1713 const auto& controls = this->well_ecl_.productionControls(summary_state);
1714 bhp = controls.bhp_limit;
1715 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1720 template<
typename TypeTag>
1722 StandardWell<TypeTag>::
1723 computeWellRatesWithThpAlqProd(
const Simulator& simulator,
1724 const SummaryState& summary_state,
1725 DeferredLogger& deferred_logger,
1726 std::vector<Scalar>& potentials,
1730 computeWellRatesAndBhpWithThpAlqProd(simulator,
1737 template<
typename TypeTag>
1742 std::vector<Scalar>& well_potentials,
1745 const auto [compute_potential, bhp_controlled_well] =
1748 if (!compute_potential) {
1752 bool converged_implicit =
false;
1753 if (this->param_.local_well_solver_control_switching_) {
1754 converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
1756 if (!converged_implicit) {
1758 const auto& summaryState = simulator.vanguard().summaryState();
1759 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1769 const auto& ws = well_state.well(this->index_of_well_);
1770 if (this->isInjector())
1771 bhp = std::max(ws.bhp, bhp);
1773 bhp = std::min(ws.bhp, bhp);
1775 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1776 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1779 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1783 this->checkNegativeWellPotentials(well_potentials,
1784 this->param_.check_well_operability_,
1794 template<
typename TypeTag>
1795 typename StandardWell<TypeTag>::Scalar
1798 const int openConnIdx)
const
1800 return (openConnIdx < 0)
1802 : this->connections_.rho(openConnIdx);
1809 template<
typename TypeTag>
1811 StandardWell<TypeTag>::
1812 updatePrimaryVariables(
const Simulator& simulator,
1813 const WellState<Scalar>& well_state,
1814 DeferredLogger& deferred_logger)
1816 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1818 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1819 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1822 if constexpr (Base::has_polymermw) {
1823 this->primary_variables_.updatePolyMW(well_state);
1826 this->primary_variables_.checkFinite(deferred_logger);
1832 template<
typename TypeTag>
1833 typename StandardWell<TypeTag>::Scalar
1834 StandardWell<TypeTag>::
1835 getRefDensity()
const
1837 return this->connections_.rho();
1843 template<
typename TypeTag>
1845 StandardWell<TypeTag>::
1846 updateWaterMobilityWithPolymer(
const Simulator& simulator,
1848 std::vector<EvalWell>& mob,
1849 DeferredLogger& deferred_logger)
const
1851 const int cell_idx = this->well_cells_[perf];
1852 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, 0);
1853 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1857 if (this->isInjector()) {
1859 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1860 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1861 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1864 if (PolymerModule::hasPlyshlog()) {
1867 if (this->isInjector() && this->wpolymer() == 0.) {
1872 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1873 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1875 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1876 PerforationRates<Scalar> perf_rates;
1877 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1878 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1879 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1880 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1881 perf_rates, deferred_logger);
1883 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1884 const auto& material_law_manager = simulator.problem().materialLawManager();
1885 const auto& scaled_drainage_info =
1886 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1887 const Scalar swcr = scaled_drainage_info.Swcr;
1888 const EvalWell poro = this->extendEval(int_quant.porosity());
1889 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1891 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1892 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1893 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1895 if (PolymerModule::hasShrate()) {
1898 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1900 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1901 int_quant.pvtRegionIndex(),
1904 mob[waterCompIdx] /= shear_factor;
1908 template<
typename TypeTag>
1910 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
1912 this->linSys_.extract(jacobian);
1916 template <
typename TypeTag>
1918 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1919 const BVector& weights,
1920 const int pressureVarIndex,
1921 const bool use_well_weights,
1922 const WellState<Scalar>& well_state)
const
1924 this->linSys_.extractCPRPressureMatrix(jacobian,
1935 template<
typename TypeTag>
1936 typename StandardWell<TypeTag>::EvalWell
1937 StandardWell<TypeTag>::
1938 pskinwater(
const Scalar throughput,
1939 const EvalWell& water_velocity,
1940 DeferredLogger& deferred_logger)
const
1942 if constexpr (Base::has_polymermw) {
1943 const int water_table_id = this->polymerWaterTable_();
1944 if (water_table_id <= 0) {
1945 OPM_DEFLOG_THROW(std::runtime_error,
1946 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
1949 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1950 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1952 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1953 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1956 OPM_DEFLOG_THROW(std::runtime_error,
1957 fmt::format(
"Polymermw is not activated, while injecting "
1958 "skin pressure is requested for well {}", name()),
1967 template<
typename TypeTag>
1968 typename StandardWell<TypeTag>::EvalWell
1969 StandardWell<TypeTag>::
1970 pskin(
const Scalar throughput,
1971 const EvalWell& water_velocity,
1972 const EvalWell& poly_inj_conc,
1973 DeferredLogger& deferred_logger)
const
1975 if constexpr (Base::has_polymermw) {
1976 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1977 const EvalWell water_velocity_abs = abs(water_velocity);
1978 if (poly_inj_conc == 0.) {
1979 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1981 const int polymer_table_id = this->polymerTable_();
1982 if (polymer_table_id <= 0) {
1983 OPM_DEFLOG_THROW(std::runtime_error,
1984 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
1987 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1988 const Scalar reference_concentration = skprpolytable.refConcentration;
1989 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1991 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1992 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1993 if (poly_inj_conc == reference_concentration) {
1994 return sign * pskin_poly;
1997 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1998 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1999 return sign * pskin;
2001 OPM_DEFLOG_THROW(std::runtime_error,
2002 fmt::format(
"Polymermw is not activated, while injecting "
2003 "skin pressure is requested for well {}", name()),
2012 template<
typename TypeTag>
2013 typename StandardWell<TypeTag>::EvalWell
2014 StandardWell<TypeTag>::
2015 wpolymermw(
const Scalar throughput,
2016 const EvalWell& water_velocity,
2017 DeferredLogger& deferred_logger)
const
2019 if constexpr (Base::has_polymermw) {
2020 const int table_id = this->polymerInjTable_();
2021 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2022 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2023 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2024 if (this->wpolymer() == 0.) {
2025 return molecular_weight;
2027 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2028 return molecular_weight;
2030 OPM_DEFLOG_THROW(std::runtime_error,
2031 fmt::format(
"Polymermw is not activated, while injecting "
2032 "polymer molecular weight is requested for well {}", name()),
2041 template<
typename TypeTag>
2043 StandardWell<TypeTag>::
2044 updateWaterThroughput([[maybe_unused]]
const double dt,
2045 WellState<Scalar>& well_state)
const
2047 if constexpr (Base::has_polymermw) {
2048 if (!this->isInjector()) {
2052 auto& perf_water_throughput = well_state.well(this->index_of_well_)
2053 .perf_data.water_throughput;
2055 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2056 const Scalar perf_water_vel =
2057 this->primary_variables_.value(Bhp + 1 + perf);
2061 if (perf_water_vel > Scalar{0}) {
2062 perf_water_throughput[perf] += perf_water_vel * dt;
2072 template<
typename TypeTag>
2074 StandardWell<TypeTag>::
2075 handleInjectivityRate(
const Simulator& simulator,
2077 std::vector<EvalWell>& cq_s)
const
2079 const int cell_idx = this->well_cells_[perf];
2080 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2081 const auto& fs = int_quants.fluidState();
2082 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2083 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2084 const int wat_vel_index = Bhp + 1 + perf;
2085 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2089 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2095 template<
typename TypeTag>
2097 StandardWell<TypeTag>::
2098 handleInjectivityEquations(
const Simulator& simulator,
2099 const WellState<Scalar>& well_state,
2101 const EvalWell& water_flux_s,
2102 DeferredLogger& deferred_logger)
2104 const int cell_idx = this->well_cells_[perf];
2105 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2106 const auto& fs = int_quants.fluidState();
2107 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2108 const EvalWell water_flux_r = water_flux_s / b_w;
2109 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2110 const EvalWell water_velocity = water_flux_r / area;
2111 const int wat_vel_index = Bhp + 1 + perf;
2114 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2116 const auto& ws = well_state.well(this->index_of_well_);
2117 const auto& perf_data = ws.perf_data;
2118 const auto& perf_water_throughput = perf_data.water_throughput;
2119 const Scalar throughput = perf_water_throughput[perf];
2120 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2122 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2123 poly_conc.setValue(this->wpolymer());
2126 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2127 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2129 StandardWellAssemble<FluidSystem,Indices>(*this).
2130 assembleInjectivityEq(eq_pskin,
2135 this->primary_variables_.numWellEq(),
2143 template<
typename TypeTag>
2145 StandardWell<TypeTag>::
2146 checkConvergenceExtraEqs(
const std::vector<Scalar>& res,
2147 ConvergenceReport& report)
const
2152 if constexpr (Base::has_polymermw) {
2153 WellConvergence(*this).
2154 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2162 template<
typename TypeTag>
2164 StandardWell<TypeTag>::
2165 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2166 const IntensiveQuantities& int_quants,
2167 const WellState<Scalar>& well_state,
2169 std::vector<RateVector>& connectionRates,
2170 DeferredLogger& deferred_logger)
const
2173 EvalWell cq_s_polymw = cq_s_poly;
2174 if (this->isInjector()) {
2175 const int wat_vel_index = Bhp + 1 + perf;
2176 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2177 if (water_velocity > 0.) {
2178 const auto& ws = well_state.well(this->index_of_well_);
2179 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2180 const Scalar throughput = perf_water_throughput[perf];
2181 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2182 cq_s_polymw *= molecular_weight;
2188 }
else if (this->isProducer()) {
2189 if (cq_s_polymw < 0.) {
2190 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2197 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2205 template<
typename TypeTag>
2206 std::optional<typename StandardWell<TypeTag>::Scalar>
2207 StandardWell<TypeTag>::
2208 computeBhpAtThpLimitProd(
const WellState<Scalar>& well_state,
2209 const Simulator& simulator,
2210 const SummaryState& summary_state,
2211 DeferredLogger& deferred_logger)
const
2213 return computeBhpAtThpLimitProdWithAlq(simulator,
2215 this->getALQ(well_state),
2219 template<
typename TypeTag>
2220 std::optional<typename StandardWell<TypeTag>::Scalar>
2221 StandardWell<TypeTag>::
2222 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2223 const SummaryState& summary_state,
2224 const Scalar alq_value,
2225 DeferredLogger& deferred_logger)
const
2228 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2234 std::vector<Scalar> rates(3);
2235 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2236 this->adaptRatesForVFP(rates);
2240 Scalar max_pressure = 0.0;
2241 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2242 const int cell_idx = this->well_cells_[perf];
2243 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2244 const auto& fs = int_quants.fluidState();
2245 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2246 max_pressure = std::max(max_pressure, pressure_cell);
2248 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2251 this->connections_.rho(),
2253 this->getTHPConstraint(summary_state),
2257 auto v = frates(*bhpAtLimit);
2258 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2264 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2268 std::vector<Scalar> rates(3);
2269 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2270 this->adaptRatesForVFP(rates);
2274 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2277 this->connections_.rho(),
2279 this->getTHPConstraint(summary_state),
2285 auto v = frates(*bhpAtLimit);
2286 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2292 return std::nullopt;
2297 template<
typename TypeTag>
2298 std::optional<typename StandardWell<TypeTag>::Scalar>
2299 StandardWell<TypeTag>::
2300 computeBhpAtThpLimitInj(
const Simulator& simulator,
2301 const SummaryState& summary_state,
2302 DeferredLogger& deferred_logger)
const
2305 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2311 std::vector<Scalar> rates(3);
2312 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2316 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2318 this->connections_.rho(),
2329 template<
typename TypeTag>
2331 StandardWell<TypeTag>::
2332 iterateWellEqWithControl(
const Simulator& simulator,
2334 const Well::InjectionControls& inj_controls,
2335 const Well::ProductionControls& prod_controls,
2336 WellState<Scalar>& well_state,
2337 const GroupState<Scalar>& group_state,
2338 DeferredLogger& deferred_logger)
2340 const int max_iter = this->param_.max_inner_iter_wells_;
2343 bool relax_convergence =
false;
2344 this->regularize_ =
false;
2346 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2348 if (it > this->param_.strict_inner_iter_wells_) {
2349 relax_convergence =
true;
2350 this->regularize_ =
true;
2353 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2355 converged = report.converged();
2361 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2368 initPrimaryVariablesEvaluation();
2369 }
while (it < max_iter);
2375 template<
typename TypeTag>
2377 StandardWell<TypeTag>::
2378 iterateWellEqWithSwitching(
const Simulator& simulator,
2380 const Well::InjectionControls& inj_controls,
2381 const Well::ProductionControls& prod_controls,
2382 WellState<Scalar>& well_state,
2383 const GroupState<Scalar>& group_state,
2384 DeferredLogger& deferred_logger,
2385 const bool fixed_control ,
2386 const bool fixed_status )
2388 const int max_iter = this->param_.max_inner_iter_wells_;
2390 bool converged =
false;
2391 bool relax_convergence =
false;
2392 this->regularize_ =
false;
2393 const auto& summary_state = simulator.vanguard().summaryState();
2398 constexpr int min_its_after_switch = 4;
2399 int its_since_last_switch = min_its_after_switch;
2400 int switch_count= 0;
2402 const auto well_status_orig = this->wellStatus_;
2403 const auto operability_orig = this->operability_status_;
2404 auto well_status_cur = well_status_orig;
2405 int status_switch_count = 0;
2407 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2408 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2410 const bool allow_switching =
2411 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2412 (!fixed_control || !fixed_status) && allow_open;
2414 bool changed =
false;
2415 bool final_check =
false;
2417 this->operability_status_.resetOperability();
2418 this->operability_status_.solvable =
true;
2420 its_since_last_switch++;
2421 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2422 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2423 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2424 inj_controls, prod_controls, wqTotal,
2425 deferred_logger, fixed_control, fixed_status);
2427 its_since_last_switch = 0;
2429 if (well_status_cur != this->wellStatus_) {
2430 well_status_cur = this->wellStatus_;
2431 status_switch_count++;
2434 if (!changed && final_check) {
2437 final_check =
false;
2441 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2443 if (it > this->param_.strict_inner_iter_wells_) {
2444 relax_convergence =
true;
2445 this->regularize_ =
true;
2448 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2450 converged = report.converged();
2454 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2456 its_since_last_switch = min_its_after_switch;
2463 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2464 initPrimaryVariablesEvaluation();
2466 }
while (it < max_iter);
2469 if (allow_switching){
2471 const bool is_stopped = this->wellIsStopped();
2472 if (this->wellHasTHPConstraints(summary_state)){
2473 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2474 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2476 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2480 this->wellStatus_ = well_status_orig;
2481 this->operability_status_ = operability_orig;
2482 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
2483 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2484 deferred_logger.debug(message);
2490 template<
typename TypeTag>
2491 std::vector<typename StandardWell<TypeTag>::Scalar>
2497 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2498 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2499 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2500 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2501 const int cell_idx = this->well_cells_[perf];
2502 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2503 std::vector<Scalar> mob(this->num_components_, 0.);
2504 getMobility(simulator, perf, mob, deferred_logger);
2505 std::vector<Scalar> cq_s(this->num_components_, 0.);
2506 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2507 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2508 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2510 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2511 cq_s, perf_rates, deferred_logger);
2512 for (
int comp = 0; comp < this->num_components_; ++comp) {
2513 well_q_s[comp] += cq_s[comp];
2516 const auto& comm = this->parallel_well_info_.communication();
2517 if (comm.size() > 1)
2519 comm.sum(well_q_s.data(), well_q_s.size());
2526 template <
typename TypeTag>
2527 std::vector<typename StandardWell<TypeTag>::Scalar>
2531 const int num_pri_vars = this->primary_variables_.numWellEq();
2532 std::vector<Scalar> retval(num_pri_vars);
2533 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2534 retval[ii] = this->primary_variables_.value(ii);
2543 template <
typename TypeTag>
2545 StandardWell<TypeTag>::
2546 setPrimaryVars(
typename std::vector<Scalar>::const_iterator it)
2548 const int num_pri_vars = this->primary_variables_.numWellEq();
2549 for (
int ii = 0; ii < num_pri_vars; ++ii) {
2550 this->primary_variables_.setValue(ii, it[ii]);
2552 return num_pri_vars;
2556 template <
typename TypeTag>
2557 typename StandardWell<TypeTag>::Eval
2558 StandardWell<TypeTag>::
2559 connectionRateEnergy(
const Scalar maxOilSaturation,
2560 const std::vector<EvalWell>& cq_s,
2561 const IntensiveQuantities& intQuants,
2562 DeferredLogger& deferred_logger)
const
2564 auto fs = intQuants.fluidState();
2566 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2567 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2572 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2573 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2574 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2575 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2576 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2579 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2580 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2585 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2587 deferred_logger.debug(
2588 fmt::format(
"Problematic d value {} obtained for well {}"
2589 " during calculateSinglePerf with rs {}"
2590 ", rv {}. Continue as if no dissolution (rs = 0) and"
2591 " vaporization (rv = 0) for this connection.",
2592 d, this->name(), fs.Rs(), fs.Rv()));
2593 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2595 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2596 cq_r_thermal = (cq_s[gasCompIdx] -
2597 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2598 (d * this->extendEval(fs.invB(phaseIdx)) );
2599 }
else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2601 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2603 (d * this->extendEval(fs.invB(phaseIdx)) );
2609 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
2611 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2612 fs.setTemperature(this->well_ecl_.inj_temperature());
2613 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
2614 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2615 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2616 paramCache.setRegionIndex(pvtRegionIdx);
2617 paramCache.setMaxOilSat(maxOilSaturation);
2618 paramCache.updatePhase(fs, phaseIdx);
2620 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2621 fs.setDensity(phaseIdx, rho);
2622 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2623 fs.setEnthalpy(phaseIdx, h);
2624 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2625 result += getValue(cq_r_thermal);
2628 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2629 result += Base::restrictEval(cq_r_thermal);
2633 return result * this->well_efficiency_factor_;
2637 template <
typename TypeTag>
2638 template<
class Value>
2640 StandardWell<TypeTag>::
2641 gasOilPerfRateInj(
const std::vector<Value>& cq_s,
2642 PerforationRates<Scalar>& perf_rates,
2645 const Value& pressure,
2647 DeferredLogger& deferred_logger)
const
2649 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2650 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2659 const Scalar d = 1.0 - getValue(rv) * getValue(rs);
2662 deferred_logger.debug(dValueError(d, this->name(),
2663 "gasOilPerfRateInj",
2668 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2671 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
2674 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2679 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
2685 template <
typename TypeTag>
2686 template<
class Value>
2688 StandardWell<TypeTag>::
2689 gasOilPerfRateProd(std::vector<Value>& cq_s,
2690 PerforationRates<Scalar>& perf_rates,
2693 const Value& rvw)
const
2695 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2696 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2697 const Value cq_sOil = cq_s[oilCompIdx];
2698 const Value cq_sGas = cq_s[gasCompIdx];
2699 const Value dis_gas = rs * cq_sOil;
2700 const Value vap_oil = rv * cq_sGas;
2702 cq_s[gasCompIdx] += dis_gas;
2703 cq_s[oilCompIdx] += vap_oil;
2706 if (this->isProducer()) {
2707 perf_rates.dis_gas = getValue(dis_gas);
2708 perf_rates.vap_oil = getValue(vap_oil);
2711 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
2712 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2713 const Value vap_wat = rvw * cq_sGas;
2714 cq_s[waterCompIdx] += vap_wat;
2715 if (this->isProducer())
2716 perf_rates.vap_wat = getValue(vap_wat);
2721 template <
typename TypeTag>
2722 template<
class Value>
2724 StandardWell<TypeTag>::
2725 gasWaterPerfRateProd(std::vector<Value>& cq_s,
2726 PerforationRates<Scalar>& perf_rates,
2728 const Value& rsw)
const
2730 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2731 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2732 const Value cq_sWat = cq_s[waterCompIdx];
2733 const Value cq_sGas = cq_s[gasCompIdx];
2734 const Value vap_wat = rvw * cq_sGas;
2735 const Value dis_gas_wat = rsw * cq_sWat;
2736 cq_s[waterCompIdx] += vap_wat;
2737 cq_s[gasCompIdx] += dis_gas_wat;
2738 if (this->isProducer()) {
2739 perf_rates.vap_wat = getValue(vap_wat);
2740 perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
2745 template <
typename TypeTag>
2746 template<
class Value>
2748 StandardWell<TypeTag>::
2749 gasWaterPerfRateInj(
const std::vector<Value>& cq_s,
2750 PerforationRates<Scalar>& perf_rates,
2753 const Value& pressure,
2754 DeferredLogger& deferred_logger)
const
2757 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2758 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2760 const Scalar dw = 1.0 - getValue(rvw) * getValue(rsw);
2763 deferred_logger.debug(dValueError(dw, this->name(),
2764 "gasWaterPerfRateInj",
2765 rsw, rvw, pressure));
2769 perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
2772 perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
2777 template <
typename TypeTag>
2778 template<
class Value>
2780 StandardWell<TypeTag>::
2781 disOilVapWatVolumeRatio(Value& volumeRatio,
2784 const Value& pressure,
2785 const std::vector<Value>& cmix_s,
2786 const std::vector<Value>& b_perfcells_dense,
2787 DeferredLogger& deferred_logger)
const
2789 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2790 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2792 const Value d = 1.0 - rvw * rsw;
2795 deferred_logger.debug(dValueError(d, this->name(),
2796 "disOilVapWatVolumeRatio",
2797 rsw, rvw, pressure));
2799 const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
2800 : cmix_s[waterCompIdx];
2801 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
2803 const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
2804 : cmix_s[gasCompIdx];
2805 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
2809 template <
typename TypeTag>
2810 template<
class Value>
2812 StandardWell<TypeTag>::
2813 gasOilVolumeRatio(Value& volumeRatio,
2816 const Value& pressure,
2817 const std::vector<Value>& cmix_s,
2818 const std::vector<Value>& b_perfcells_dense,
2819 DeferredLogger& deferred_logger)
const
2821 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2822 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2824 const Value d = 1.0 - rv * rs;
2827 deferred_logger.debug(dValueError(d, this->name(),
2828 "gasOilVolumeRatio",
2831 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
2832 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
2834 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
2835 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
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:291
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:278
Definition StandardWell.hpp:59
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition StandardWell_impl.hpp:1417
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 WellTest.hpp:37
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:62
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
Definition PerforationData.hpp:41