22 #include <opm/common/utility/numeric/RootFinders.hpp>
23 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellInjectionProperties.hpp>
24 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
25 #include <opm/simulators/linalg/MatrixBlock.hpp>
26 #include <opm/simulators/wells/VFPHelpers.hpp>
35 template<
typename TypeTag>
36 StandardWell<TypeTag>::
37 StandardWell(
const Well& well,
38 const ParallelWellInfo& pw_info,
40 const ModelParameters& param,
41 const RateConverterType& rate_converter,
42 const int pvtRegionIdx,
43 const int num_components,
45 const int index_of_well,
46 const std::vector<PerforationData>& perf_data)
47 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
48 , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
50 assert(this->num_components_ == numWellConservationEq);
57 template<
typename TypeTag>
59 StandardWell<TypeTag>::
60 init(
const PhaseUsage* phase_usage_arg,
61 const std::vector<double>& depth_arg,
62 const double gravity_arg,
64 const std::vector< Scalar >& B_avg)
66 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
67 this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
74 template<
typename TypeTag>
75 void StandardWell<TypeTag>::
76 initPrimaryVariablesEvaluation()
const
78 this->StdWellEval::initPrimaryVariablesEvaluation();
85 template<
typename TypeTag>
86 typename StandardWell<TypeTag>::Eval
87 StandardWell<TypeTag>::getPerfCellPressure(
const typename StandardWell<TypeTag>::FluidState& fs)
const
90 if (Indices::oilEnabled) {
91 pressure = fs.pressure(FluidSystem::oilPhaseIdx);
93 if (Indices::waterEnabled) {
94 pressure = fs.pressure(FluidSystem::waterPhaseIdx);
96 pressure = fs.pressure(FluidSystem::gasPhaseIdx);
103 template<
typename TypeTag>
105 StandardWell<TypeTag>::
106 computePerfRateEval(
const IntensiveQuantities& intQuants,
107 const std::vector<EvalWell>& mob,
112 std::vector<EvalWell>& cq_s,
113 double& perf_dis_gas_rate,
114 double& perf_vap_oil_rate,
115 DeferredLogger& deferred_logger)
const
117 const auto& fs = intQuants.fluidState();
118 const EvalWell pressure = this->extendEval(getPerfCellPressure(fs));
119 const EvalWell rs = this->extendEval(fs.Rs());
120 const EvalWell rv = this->extendEval(fs.Rv());
121 std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
122 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
123 if (!FluidSystem::phaseIsActive(phaseIdx)) {
126 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
127 b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
129 if constexpr (has_solvent) {
130 b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
133 if constexpr (has_zFraction) {
134 if (this->isInjector()) {
135 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
136 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
137 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
141 EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
143 if (this->isInjector()) {
144 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
145 skin_pressure = this->primary_variables_evaluation_[pskin_index];
150 std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
151 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
152 cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
172 template<
typename TypeTag>
174 StandardWell<TypeTag>::
175 computePerfRateScalar(
const IntensiveQuantities& intQuants,
176 const std::vector<Scalar>& mob,
181 std::vector<Scalar>& cq_s,
182 DeferredLogger& deferred_logger)
const
184 const auto& fs = intQuants.fluidState();
185 const Scalar pressure = getPerfCellPressure(fs).value();
186 const Scalar rs = fs.Rs().value();
187 const Scalar rv = fs.Rv().value();
188 std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
189 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
190 if (!FluidSystem::phaseIsActive(phaseIdx)) {
193 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
194 b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
196 if constexpr (has_solvent) {
197 b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
200 if constexpr (has_zFraction) {
201 if (this->isInjector()) {
202 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
203 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
204 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
208 Scalar skin_pressure =0.0;
210 if (this->isInjector()) {
211 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
212 skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
216 Scalar perf_dis_gas_rate = 0.0;
217 Scalar perf_vap_oil_rate = 0.0;
220 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
221 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
222 cmix_s[componentIdx] = getValue(this->wellSurfaceVolumeFraction(componentIdx));
242 template<
typename TypeTag>
243 template<
class Value>
245 StandardWell<TypeTag>::
246 computePerfRate(
const std::vector<Value>& mob,
247 const Value& pressure,
251 std::vector<Value>& b_perfcells_dense,
255 const Value& skin_pressure,
256 const std::vector<Value>& cmix_s,
257 std::vector<Value>& cq_s,
258 double& perf_dis_gas_rate,
259 double& perf_vap_oil_rate,
260 DeferredLogger& deferred_logger)
const
263 const Value well_pressure = bhp + this->perf_pressure_diffs_[perf];
264 Value drawdown = pressure - well_pressure;
265 if (this->isInjector()) {
266 drawdown += skin_pressure;
270 if ( drawdown > 0 ) {
272 if (!allow_cf && this->isInjector()) {
277 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
278 const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
279 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
282 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
283 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
284 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
285 const Value cq_sOil = cq_s[oilCompIdx];
286 const Value cq_sGas = cq_s[gasCompIdx];
287 const Value dis_gas = rs * cq_sOil;
288 const Value vap_oil = rv * cq_sGas;
290 cq_s[gasCompIdx] += dis_gas;
291 cq_s[oilCompIdx] += vap_oil;
294 if (this->isProducer()) {
295 perf_dis_gas_rate = getValue(dis_gas);
296 perf_vap_oil_rate = getValue(vap_oil);
302 if (!allow_cf && this->isProducer()) {
307 Value total_mob_dense = mob[0];
308 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
309 total_mob_dense += mob[componentIdx];
313 const Value cqt_i = - Tw * (total_mob_dense * drawdown);
316 Value volumeRatio = bhp * 0.0;
318 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
319 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
320 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
323 if constexpr (Indices::enableSolvent) {
324 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
327 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
328 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
329 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
331 const Value d = 1.0 - rv * rs;
333 if (getValue(d) == 0.0) {
334 OPM_DEFLOG_THROW(NumericalIssue,
"Zero d value obtained for well " << this->name() <<
" during flux calcuation"
335 <<
" with rs " << rs <<
" and rv " << rv, deferred_logger);
338 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
339 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
341 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
342 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
345 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
346 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
347 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
349 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
350 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
351 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
356 Value cqt_is = cqt_i/volumeRatio;
357 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
358 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
362 if (this->isProducer()) {
363 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
364 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
365 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
374 const double d = 1.0 - getValue(rv) * getValue(rs);
377 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
380 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
387 template<
typename TypeTag>
389 StandardWell<TypeTag>::
390 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
392 const Well::InjectionControls& ,
393 const Well::ProductionControls& ,
394 WellState& well_state,
395 const GroupState& group_state,
396 DeferredLogger& deferred_logger)
400 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
405 this->invDuneD_ = 0.0;
406 this->resWell_ = 0.0;
408 assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
414 template<
typename TypeTag>
416 StandardWell<TypeTag>::
417 assembleWellEqWithoutIterationImpl(
const Simulator& ebosSimulator,
419 WellState& well_state,
420 const GroupState& group_state,
421 DeferredLogger& deferred_logger)
425 const double volume = 0.002831684659200;
426 auto& ws = well_state.well(this->index_of_well_);
428 ws.vaporized_oil_rate = 0;
429 ws.dissolved_gas_rate = 0;
431 const int np = this->number_of_phases_;
433 std::vector<RateVector> connectionRates = this->connectionRates_;
434 auto& perf_data = ws.perf_data;
435 auto& perf_rates = perf_data.phase_rates;
436 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
438 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
439 EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
440 EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
441 calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
444 if constexpr (has_polymer && Base::has_polymermw) {
445 if (this->isInjector()) {
446 handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
449 const int cell_idx = this->well_cells_[perf];
450 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
452 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
454 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
457 this->resWell_[0][componentIdx] += cq_s_effective.value();
460 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
462 this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq);
463 this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
466 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
467 this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
471 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
472 auto& perf_rate_solvent = perf_data.solvent_rates;
473 perf_rate_solvent[perf] = cq_s[componentIdx].value();
475 perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
479 if constexpr (has_zFraction) {
480 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
481 this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
486 this->connectionRates_ = connectionRates;
489 wellhelpers::sumDistributedWellEntries(this->invDuneD_[0][0], this->resWell_[0],
490 this->parallel_well_info_.communication());
492 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
495 EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
496 if (FluidSystem::numActivePhases() > 1) {
498 resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
500 resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
501 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
502 this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
504 this->resWell_[0][componentIdx] += resWell_loc.value();
507 const auto& summaryState = ebosSimulator.vanguard().summaryState();
508 const Schedule& schedule = ebosSimulator.vanguard().schedule();
509 this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger);
514 Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]);
516 OPM_DEFLOG_THROW(NumericalIssue,
"Error when inverting local well equations for well " + name(), deferred_logger);
523 template<
typename TypeTag>
525 StandardWell<TypeTag>::
526 calculateSinglePerf(
const Simulator& ebosSimulator,
528 WellState& well_state,
529 std::vector<RateVector>& connectionRates,
530 std::vector<EvalWell>& cq_s,
531 EvalWell& water_flux_s,
532 EvalWell& cq_s_zfrac_effective,
533 DeferredLogger& deferred_logger)
const
535 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
536 const EvalWell& bhp = this->getBhp();
537 const int cell_idx = this->well_cells_[perf];
538 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
539 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
540 getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
542 double perf_dis_gas_rate = 0.;
543 double perf_vap_oil_rate = 0.;
544 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
545 const double Tw = this->well_index_[perf] * trans_mult;
546 computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
547 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
549 auto& ws = well_state.well(this->index_of_well_);
550 auto& perf_data = ws.perf_data;
551 if constexpr (has_polymer && Base::has_polymermw) {
552 if (this->isInjector()) {
555 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
556 water_flux_s = cq_s[water_comp_idx];
559 handleInjectivityRate(ebosSimulator, perf, cq_s);
564 if (this->isProducer()) {
565 ws.dissolved_gas_rate += perf_dis_gas_rate;
566 ws.vaporized_oil_rate += perf_vap_oil_rate;
569 if constexpr (has_energy) {
570 connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
573 if constexpr (has_energy) {
575 auto fs = intQuants.fluidState();
576 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
577 if (!FluidSystem::phaseIsActive(phaseIdx)) {
581 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
583 EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
584 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
586 if(FluidSystem::waterPhaseIdx == phaseIdx)
587 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
590 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
591 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
595 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
597 if(FluidSystem::gasPhaseIdx == phaseIdx)
598 cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
600 if(FluidSystem::oilPhaseIdx == phaseIdx)
601 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
604 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
608 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
610 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
611 fs.setTemperature(this->well_ecl_.temperature());
612 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
613 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
614 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
615 paramCache.setRegionIndex(pvtRegionIdx);
616 paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
617 paramCache.updatePhase(fs, phaseIdx);
619 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
620 fs.setDensity(phaseIdx, rho);
621 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
622 fs.setEnthalpy(phaseIdx, h);
625 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
626 connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
630 if constexpr (has_polymer) {
632 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
633 EvalWell cq_s_poly = cq_s[waterCompIdx];
634 if (this->isInjector()) {
635 cq_s_poly *= this->wpolymer();
637 cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
640 auto& perf_rate_polymer = perf_data.polymer_rates;
641 perf_rate_polymer[perf] = cq_s_poly.value();
643 cq_s_poly *= this->well_efficiency_factor_;
644 connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
646 if constexpr (Base::has_polymermw) {
647 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
651 if constexpr (has_foam) {
653 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
654 EvalWell cq_s_foam = cq_s[gasCompIdx] * this->well_efficiency_factor_;
655 if (this->isInjector()) {
656 cq_s_foam *= this->wfoam();
658 cq_s_foam *= this->extendEval(intQuants.foamConcentration());
660 connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
663 if constexpr (has_zFraction) {
665 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
666 cq_s_zfrac_effective = cq_s[gasCompIdx];
667 if (this->isInjector()) {
668 cq_s_zfrac_effective *= this->wsolvent();
669 }
else if (cq_s_zfrac_effective.value() != 0.0) {
670 const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
671 cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
673 auto& perf_rate_solvent = perf_data.solvent_rates;
674 perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
676 cq_s_zfrac_effective *= this->well_efficiency_factor_;
677 connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
680 if constexpr (has_brine) {
682 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
683 EvalWell cq_s_sm = cq_s[waterCompIdx];
684 if (this->isInjector()) {
685 cq_s_sm *= this->wsalt();
687 cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration());
690 auto& perf_rate_brine = perf_data.brine_rates;
691 perf_rate_brine[perf] = cq_s_sm.value();
693 cq_s_sm *= this->well_efficiency_factor_;
694 connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
697 if constexpr (has_micp) {
698 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
699 EvalWell cq_s_microbe = cq_s[waterCompIdx];
700 if (this->isInjector()) {
701 cq_s_microbe *= this->wmicrobes();
703 cq_s_microbe *= this->extendEval(intQuants.microbialConcentration());
705 connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe);
706 EvalWell cq_s_oxygen = cq_s[waterCompIdx];
707 if (this->isInjector()) {
708 cq_s_oxygen *= this->woxygen();
710 cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration());
712 connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen);
713 EvalWell cq_s_urea = cq_s[waterCompIdx];
714 if (this->isInjector()) {
715 cq_s_urea *= this->wurea();
717 cq_s_urea *= this->extendEval(intQuants.ureaConcentration());
719 connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea);
723 perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf];
729 template<
typename TypeTag>
731 StandardWell<TypeTag>::
732 getMobilityEval(
const Simulator& ebosSimulator,
734 std::vector<EvalWell>& mob,
735 DeferredLogger& deferred_logger)
const
737 const int cell_idx = this->well_cells_[perf];
738 assert (
int(mob.size()) == this->num_components_);
739 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
740 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
744 const int satid = this->saturation_table_number_[perf] - 1;
745 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
746 if( satid == satid_elem ) {
748 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
749 if (!FluidSystem::phaseIsActive(phaseIdx)) {
753 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
754 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
757 mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
761 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
762 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
763 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
766 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
769 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
770 if (!FluidSystem::phaseIsActive(phaseIdx)) {
774 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
775 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
779 if constexpr (has_solvent) {
780 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
785 if constexpr (has_polymer) {
786 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
787 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
792 if constexpr (!Base::has_polymermw) {
793 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
798 template<
typename TypeTag>
800 StandardWell<TypeTag>::
801 getMobilityScalar(
const Simulator& ebosSimulator,
803 std::vector<Scalar>& mob,
804 DeferredLogger& deferred_logger)
const
806 const int cell_idx = this->well_cells_[perf];
807 assert (
int(mob.size()) == this->num_components_);
808 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
809 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
813 const int satid = this->saturation_table_number_[perf] - 1;
814 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
815 if( satid == satid_elem ) {
817 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
818 if (!FluidSystem::phaseIsActive(phaseIdx)) {
822 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
823 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
826 mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
830 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
831 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
832 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
835 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
838 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
839 if (!FluidSystem::phaseIsActive(phaseIdx)) {
843 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
844 mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx));
848 if constexpr (has_solvent) {
849 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
854 if constexpr (has_polymer) {
855 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
856 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
861 if constexpr (!Base::has_polymermw) {
862 std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
863 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
864 for (
size_t i = 0; i < mob.size(); ++i) {
865 mob[i] = getValue(mob_eval[i]);
873 template<
typename TypeTag>
875 StandardWell<TypeTag>::
876 updateWellState(
const BVectorWell& dwells,
877 WellState& well_state,
878 DeferredLogger& deferred_logger)
const
880 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
882 updatePrimaryVariablesNewton(dwells, well_state, deferred_logger);
884 updateWellStateFromPrimaryVariables(well_state, deferred_logger);
885 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
892 template<
typename TypeTag>
894 StandardWell<TypeTag>::
895 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
897 DeferredLogger& deferred_logger)
const
899 const double dFLimit = this->param_.dwell_fraction_max_;
900 const double dBHPLimit = this->param_.dbhp_max_rel_;
901 this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
903 updateExtraPrimaryVariables(dwells);
905 for (
double v : this->primary_variables_) {
907 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after newton update well: " << this->name(), deferred_logger);
916 template<
typename TypeTag>
918 StandardWell<TypeTag>::
919 updateExtraPrimaryVariables(
const BVectorWell& dwells)
const
922 if constexpr (Base::has_polymermw) {
923 this->updatePrimaryVariablesPolyMW(dwells);
931 template<
typename TypeTag>
933 StandardWell<TypeTag>::
934 updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger)
const
936 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger);
939 if constexpr (Base::has_polymermw) {
940 this->updateWellStateFromPrimaryVariablesPolyMW(well_state);
948 template<
typename TypeTag>
950 StandardWell<TypeTag>::
951 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
957 if (this->isInjector()) {
962 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
963 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
965 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
966 std::vector<Scalar> mob(this->num_components_, 0.0);
967 getMobilityScalar(ebos_simulator, perf, mob, deferred_logger);
969 const int cell_idx = this->well_cells_[perf];
970 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
971 const auto& fs = int_quantities.fluidState();
973 Eval perf_pressure = getPerfCellPressure(fs);
974 double p_r = perf_pressure.value();
977 std::vector<double> b_perf(this->num_components_);
978 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
979 if (!FluidSystem::phaseIsActive(phase)) {
982 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
983 b_perf[comp_idx] = fs.invB(phase).value();
987 const double h_perf = this->perf_pressure_diffs_[perf];
988 const double pressure_diff = p_r - h_perf;
993 if (pressure_diff <= 0.) {
994 deferred_logger.warning(
"NON_POSITIVE_DRAWDOWN_IPR",
995 "non-positive drawdown found when updateIPR for well " + name());
999 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1004 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1005 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1006 for (
int p = 0; p < this->number_of_phases_; ++p) {
1007 const double tw_mob = tw_perf * mob[p] * b_perf[p];
1008 ipr_a_perf[p] += tw_mob * pressure_diff;
1009 ipr_b_perf[p] += tw_mob;
1013 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1014 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1015 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1016 const double rs = (fs.Rs()).value();
1017 const double rv = (fs.Rv()).value();
1019 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1020 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1022 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1023 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1025 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1026 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1028 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1029 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1032 for (
int p = 0; p < this->number_of_phases_; ++p) {
1034 this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
1035 this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
1038 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1039 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1043 template<
typename TypeTag>
1045 StandardWell<TypeTag>::
1046 checkOperabilityUnderBHPLimitProducer(
const WellState& well_state,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1048 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1049 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1052 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1053 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1057 for (
int p = 0; p < this->number_of_phases_; ++p) {
1058 const double temp = this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
1060 this->operability_status_.operable_under_only_bhp_limit =
false;
1066 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1070 std::vector<double> well_rates_bhp_limit;
1071 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1073 const double thp = this->calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
1074 const double thp_limit = this->getTHPConstraint(summaryState);
1076 if (thp < thp_limit) {
1077 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1088 this->operability_status_.operable_under_only_bhp_limit =
true;
1089 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1097 template<
typename TypeTag>
1099 StandardWell<TypeTag>::
1100 checkOperabilityUnderTHPLimitProducer(
const Simulator& ebos_simulator,
const WellState& well_state, DeferredLogger& deferred_logger)
1102 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1103 const auto obtain_bhp = computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger);
1106 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1108 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1109 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1111 const double thp_limit = this->getTHPConstraint(summaryState);
1112 if (*obtain_bhp < thp_limit) {
1113 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1114 +
" bars is SMALLER than thp limit "
1115 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1116 +
" bars as a producer for well " + name();
1117 deferred_logger.debug(msg);
1120 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1121 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1122 if (!this->wellIsStopped()) {
1123 const double thp_limit = this->getTHPConstraint(summaryState);
1124 deferred_logger.debug(
" could not find bhp value at thp limit "
1125 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1126 +
" bar for well " + name() +
", the well might need to be closed ");
1135 template<
typename TypeTag>
1137 StandardWell<TypeTag>::
1138 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1140 bool all_drawdown_wrong_direction =
true;
1142 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1143 const int cell_idx = this->well_cells_[perf];
1144 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1145 const auto& fs = intQuants.fluidState();
1147 const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
1148 const double bhp = this->getBhp().value();
1151 const double well_pressure = bhp + this->perf_pressure_diffs_[perf];
1152 const double drawdown = pressure - well_pressure;
1157 if ( (drawdown < 0. && this->isInjector()) ||
1158 (drawdown > 0. && this->isProducer()) ) {
1159 all_drawdown_wrong_direction =
false;
1164 const auto& comm = this->parallel_well_info_.communication();
1165 if (comm.size() > 1)
1167 all_drawdown_wrong_direction =
1168 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1171 return all_drawdown_wrong_direction;
1177 template<
typename TypeTag>
1179 StandardWell<TypeTag>::
1180 canProduceInjectWithCurrentBhp(
const Simulator& ebos_simulator,
1181 const WellState& well_state,
1182 DeferredLogger& deferred_logger)
1184 const double bhp = well_state.well(this->index_of_well_).bhp;
1185 std::vector<double> well_rates;
1186 computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
1188 const double sign = (this->isProducer()) ? -1. : 1.;
1189 const double threshold = sign * std::numeric_limits<double>::min();
1191 bool can_produce_inject =
false;
1192 for (
const auto value : well_rates) {
1193 if (this->isProducer() && value < threshold) {
1194 can_produce_inject =
true;
1196 }
else if (this->isInjector() && value > threshold) {
1197 can_produce_inject =
true;
1202 if (!can_produce_inject) {
1203 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1206 return can_produce_inject;
1213 template<
typename TypeTag>
1215 StandardWell<TypeTag>::
1216 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1218 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1224 template<
typename TypeTag>
1226 StandardWell<TypeTag>::
1227 computePropertiesForWellConnectionPressures(
const Simulator& ebosSimulator,
1228 const WellState& well_state,
1229 std::vector<double>& b_perf,
1230 std::vector<double>& rsmax_perf,
1231 std::vector<double>& rvmax_perf,
1232 std::vector<double>& surf_dens_perf)
const
1234 const int nperf = this->number_of_perforations_;
1236 b_perf.resize(nperf * this->num_components_);
1237 surf_dens_perf.resize(nperf * this->num_components_);
1238 const auto& ws = well_state.well(this->index_of_well_);
1240 const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1241 const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1242 const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1245 if (oilPresent && gasPresent) {
1246 rsmax_perf.resize(nperf);
1247 rvmax_perf.resize(nperf);
1251 const auto& perf_press = ws.perf_data.pressure;
1252 auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp,
1256 for (
int perf = 0; perf < nperf; ++perf) {
1257 const int cell_idx = this->well_cells_[perf];
1258 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1259 const auto& fs = intQuants.fluidState();
1261 const double p_avg = (perf_press[perf] + p_above[perf])/2;
1262 const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1263 const double saltConcentration = fs.saltConcentration().value();
1266 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1267 b_perf[ waterCompIdx + perf * this->num_components_] =
1268 FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration);
1272 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1273 const int gaspos = gasCompIdx + perf * this->num_components_;
1276 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1277 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1279 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1282 rv = oilrate / gasrate;
1284 rv = std::min(rv, rvmax_perf[perf]);
1286 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
1289 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1293 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1298 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1299 const int oilpos = oilCompIdx + perf * this->num_components_;
1301 rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1302 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1304 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1307 rs = gasrate / oilrate;
1309 rs = std::min(rs, rsmax_perf[perf]);
1310 b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1312 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1315 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1320 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1321 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1325 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1326 surf_dens_perf[this->num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() );
1330 if constexpr (has_solvent) {
1331 b_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1332 surf_dens_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity();
1341 template<
typename TypeTag>
1345 const std::vector<double>& B_avg,
1347 const bool relax_tolerance)
const
1351 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1353 std::vector<double> res;
1356 this->param_.max_residual_allowed_,
1357 this->param_.tolerance_wells_,
1358 this->param_.relaxed_tolerance_flow_well_,
1362 checkConvergenceExtraEqs(res, report);
1371 template<
typename TypeTag>
1379 auto fluidState = [&ebosSimulator,
this](
const int perf)
1381 const auto cell_idx = this->well_cells_[perf];
1382 return ebosSimulator.model()
1383 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
1386 const int np = this->number_of_phases_;
1387 auto setToZero = [np](
double* x) ->
void
1389 std::fill_n(x, np, 0.0);
1392 auto addVector = [np](
const double* src,
double* dest) ->
void
1394 std::transform(src, src + np, dest, dest, std::plus<>{});
1397 auto& ws = well_state.well(this->index_of_well_);
1398 auto& perf_data = ws.perf_data;
1399 auto* wellPI = ws.productivity_index.data();
1400 auto* connPI = perf_data.prod_index.data();
1404 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1405 auto subsetPerfID = 0;
1407 for (
const auto& perf : *this->perf_data_) {
1408 auto allPerfID = perf.ecl_index;
1410 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1415 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
1416 getMobilityEval(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1418 const auto& fs = fluidState(subsetPerfID);
1421 if (this->isInjector()) {
1422 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1423 mob, connPI, deferred_logger);
1426 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1429 addVector(connPI, wellPI);
1436 const auto& comm = this->parallel_well_info_.communication();
1437 if (comm.size() > 1) {
1438 comm.sum(wellPI, np);
1441 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1442 "Internal logic error in processing connections for PI/II");
1447 template<
typename TypeTag>
1449 StandardWell<TypeTag>::
1450 computeWellConnectionDensitesPressures(
const Simulator& ebosSimulator,
1451 const WellState& well_state,
1452 const std::vector<double>& b_perf,
1453 const std::vector<double>& rsmax_perf,
1454 const std::vector<double>& rvmax_perf,
1455 const std::vector<double>& surf_dens_perf)
1458 const int nperf = this->number_of_perforations_;
1459 const int np = this->number_of_phases_;
1460 std::vector<double> perfRates(b_perf.size(),0.0);
1461 const auto& ws = well_state.well(this->index_of_well_);
1462 const auto& perf_data = ws.perf_data;
1463 const auto& perf_rates_state = perf_data.phase_rates;
1465 for (
int perf = 0; perf < nperf; ++perf) {
1466 for (
int comp = 0; comp < np; ++comp) {
1467 perfRates[perf * this->num_components_ + comp] = perf_rates_state[perf * np + this->ebosCompIdxToFlowCompIdx(comp)];
1471 if constexpr (has_solvent) {
1472 const auto& solvent_perf_rates_state = perf_data.solvent_rates;
1473 for (
int perf = 0; perf < nperf; ++perf) {
1474 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = solvent_perf_rates_state[perf];
1481 bool all_zero = std::all_of(perfRates.begin(), perfRates.end(), [](
double val) { return val == 0.0; });
1482 if ( all_zero && this->isProducer() ) {
1483 double total_tw = 0;
1484 for (
int perf = 0; perf < nperf; ++perf) {
1485 total_tw += this->well_index_[perf];
1487 for (
int perf = 0; perf < nperf; ++perf) {
1488 const int cell_idx = this->well_cells_[perf];
1489 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1490 const auto& fs = intQuants.fluidState();
1491 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1492 double total_mobility = 0.0;
1493 for (
int p = 0; p < np; ++p) {
1494 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1495 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1497 if constexpr (has_solvent) {
1498 total_mobility += intQuants.solventInverseFormationVolumeFactor().value() * intQuants.solventMobility().value();
1500 for (
int p = 0; p < np; ++p) {
1501 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1502 perfRates[perf * this->num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1504 if constexpr (has_solvent) {
1505 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility;
1510 this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1512 this->computeConnectionPressureDelta();
1519 template<
typename TypeTag>
1521 StandardWell<TypeTag>::
1522 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1523 const WellState& well_state)
1528 std::vector<double> b_perf;
1529 std::vector<double> rsmax_perf;
1530 std::vector<double> rvmax_perf;
1531 std::vector<double> surf_dens_perf;
1532 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1533 computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1540 template<
typename TypeTag>
1542 StandardWell<TypeTag>::
1543 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
1545 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1549 BVectorWell dx_well(1);
1550 dx_well[0].resize(this->numWellEq_);
1551 this->invDuneD_.mv(this->resWell_, dx_well);
1553 updateWellState(dx_well, well_state, deferred_logger);
1560 template<
typename TypeTag>
1562 StandardWell<TypeTag>::
1563 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1564 const WellState& well_state,
1565 DeferredLogger& deferred_logger)
1567 updatePrimaryVariables(well_state, deferred_logger);
1568 initPrimaryVariablesEvaluation();
1569 computeWellConnectionPressures(ebosSimulator, well_state);
1570 this->computeAccumWell();
1575 template<
typename TypeTag>
1578 apply(
const BVector& x, BVector& Ax)
const
1580 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1582 if (this->param_.matrix_add_well_contributions_)
1587 assert( this->Bx_.size() == this->duneB_.N() );
1588 assert( this->invDrw_.size() == this->invDuneD_.N() );
1591 this->parallelB_.mv(x, this->Bx_);
1596 BVectorWell& invDBx = this->invDrw_;
1597 this->invDuneD_.mv(this->Bx_, invDBx);
1600 this->duneC_.mmtv(invDBx,Ax);
1606 template<
typename TypeTag>
1609 apply(BVector& r)
const
1611 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1613 assert( this->invDrw_.size() == this->invDuneD_.N() );
1616 this->invDuneD_.mv(this->resWell_, this->invDrw_);
1618 this->duneC_.mmtv(this->invDrw_, r);
1621 template<
typename TypeTag>
1626 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1628 BVectorWell resWell = this->resWell_;
1630 this->parallelB_.mmv(x, resWell);
1632 this->invDuneD_.mv(resWell, xw);
1639 template<
typename TypeTag>
1646 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1649 xw[0].resize(this->numWellEq_);
1651 recoverSolutionWell(x, xw);
1652 updateWellState(xw, well_state, deferred_logger);
1658 template<
typename TypeTag>
1663 std::vector<double>& well_flux,
1667 const int np = this->number_of_phases_;
1668 well_flux.resize(np, 0.0);
1670 const bool allow_cf = this->getAllowCrossFlow();
1672 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1673 const int cell_idx = this->well_cells_[perf];
1674 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1676 std::vector<Scalar> mob(this->num_components_, 0.);
1677 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
1678 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
1679 const double Tw = this->well_index_[perf] * trans_mult;
1681 std::vector<Scalar> cq_s(this->num_components_, 0.);
1682 computePerfRateScalar(intQuants, mob, bhp, Tw, perf, allow_cf,
1683 cq_s, deferred_logger);
1685 for(
int p = 0; p < np; ++p) {
1686 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
1689 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1694 template<
typename TypeTag>
1696 StandardWell<TypeTag>::
1697 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
1699 std::vector<double>& well_flux,
1700 DeferredLogger& deferred_logger)
const
1706 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1707 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1708 auto& ws = well_state_copy.well(this->index_of_well_);
1711 if (this->well_ecl_.isInjector()) {
1712 ws.injection_cmode = Well::InjectorCMode::BHP;
1714 ws.production_cmode = Well::ProducerCMode::BHP;
1719 const int np = this->number_of_phases_;
1720 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1721 for (
int phase = 0; phase < np; ++phase){
1722 well_state_copy.wellRates(this->index_of_well_)[phase]
1723 = sign * ws.well_potentials[phase];
1727 StandardWell<TypeTag> well(*
this);
1728 well.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1730 const double dt = ebosSimulator.timeStepSize();
1731 bool converged = well.iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
1733 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1734 " potentials are computed based on unconverged solution";
1735 deferred_logger.debug(msg);
1737 well.updatePrimaryVariables(well_state_copy, deferred_logger);
1738 well.computeWellConnectionPressures(ebosSimulator, well_state_copy);
1739 well.initPrimaryVariablesEvaluation();
1740 well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1746 template<
typename TypeTag>
1748 StandardWell<TypeTag>::
1749 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
1750 DeferredLogger& deferred_logger,
1751 const WellState &well_state)
const
1753 std::vector<double> potentials(this->number_of_phases_, 0.0);
1754 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1756 const auto& well = this->well_ecl_;
1757 if (well.isInjector()){
1758 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1759 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1760 if (bhp_at_thp_limit) {
1761 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1762 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
1764 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1765 "Failed in getting converged thp based potential calculation for well "
1766 + name() +
". Instead the bhp based value is used");
1767 const double bhp = controls.bhp_limit;
1768 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
1771 computeWellRatesWithThpAlqProd(
1772 ebos_simulator, summary_state,
1773 deferred_logger, potentials, this->getALQ(well_state)
1780 template<
typename TypeTag>
1782 StandardWell<TypeTag>::
1783 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &ebos_simulator,
1784 const SummaryState &summary_state,
1785 DeferredLogger &deferred_logger,
1786 std::vector<double> &potentials,
1790 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1791 ebos_simulator, summary_state, deferred_logger, alq);
1792 if (bhp_at_thp_limit) {
1793 const auto& controls = this->well_ecl_.productionControls(summary_state);
1794 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1795 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
1798 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1799 "Failed in getting converged thp based potential calculation for well "
1800 + name() +
". Instead the bhp based value is used");
1801 const auto& controls = this->well_ecl_.productionControls(summary_state);
1802 bhp = controls.bhp_limit;
1803 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
1808 template<
typename TypeTag>
1810 StandardWell<TypeTag>::
1811 computeWellRatesWithThpAlqProd(
const Simulator &ebos_simulator,
1812 const SummaryState &summary_state,
1813 DeferredLogger &deferred_logger,
1814 std::vector<double> &potentials,
1818 computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1825 template<
typename TypeTag>
1827 StandardWell<TypeTag>::
1828 gasLiftOptimizationStage1(
1829 WellState& well_state,
1830 const GroupState& group_state,
1831 const Simulator& ebos_simulator,
1832 DeferredLogger& deferred_logger,
1833 GLiftProdWells &prod_wells,
1834 GLiftOptWells &glift_wells,
1835 GLiftWellStateMap &glift_state_map,
1836 GasLiftGroupInfo &group_info,
1837 GLiftSyncGroups &sync_groups
1840 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1841 std::unique_ptr<GasLiftSingleWell> glift
1842 = std::make_unique<GasLiftSingleWell>(
1843 *
this, ebos_simulator, summary_state,
1844 deferred_logger, well_state, group_state, group_info, sync_groups);
1845 auto state = glift->runOptimize(
1846 ebos_simulator.model().newtonMethod().numIterations());
1848 glift_state_map.insert({this->name(), std::move(state)});
1849 glift_wells.insert({this->name(), std::move(glift)});
1852 prod_wells.insert({this->name(),
this});
1855 template<
typename TypeTag>
1860 std::vector<double>& well_potentials,
1863 const int np = this->number_of_phases_;
1864 well_potentials.resize(np, 0.0);
1866 if (this->wellIsStopped()) {
1871 bool thp_controlled_well =
false;
1872 bool bhp_controlled_well =
false;
1873 const auto& ws = well_state.well(this->index_of_well_);
1874 if (this->isInjector()) {
1875 const Well::InjectorCMode& current = ws.injection_cmode;
1876 if (current == Well::InjectorCMode::THP) {
1877 thp_controlled_well =
true;
1879 if (current == Well::InjectorCMode::BHP) {
1880 bhp_controlled_well =
true;
1883 const Well::ProducerCMode& current = ws.production_cmode;
1884 if (current == Well::ProducerCMode::THP) {
1885 thp_controlled_well =
true;
1887 if (current == Well::ProducerCMode::BHP) {
1888 bhp_controlled_well =
true;
1891 if (thp_controlled_well || bhp_controlled_well) {
1893 double total_rate = 0.0;
1894 for (
int phase = 0; phase < np; ++phase){
1895 total_rate += ws.surface_rates[phase];
1900 if (std::abs(total_rate) > 0) {
1901 for (
int phase = 0; phase < np; ++phase){
1902 well_potentials[phase] = ws.surface_rates[phase];
1909 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1910 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1912 const double bhp = this->mostStrictBhpFromBhpLimits(summaryState);
1913 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1914 computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
1917 well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
1925 template<
typename TypeTag>
1930 this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
1931 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1934 if constexpr (Base::has_polymermw) {
1935 if (this->isInjector()) {
1936 const auto& ws = well_state.well(this->index_of_well_);
1937 const auto& perf_data = ws.perf_data;
1938 const auto& water_velocity = perf_data.water_velocity;
1939 const auto& skin_pressure = perf_data.skin_pressure;
1940 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1941 this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
1942 this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
1946 for (
double v : this->primary_variables_) {
1948 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after update from wellState well: " << this->name(), deferred_logger);
1955 template<
typename TypeTag>
1957 StandardWell<TypeTag>::
1958 getRefDensity()
const
1960 return this->perf_densities_[0];
1966 template<
typename TypeTag>
1968 StandardWell<TypeTag>::
1969 updateWaterMobilityWithPolymer(
const Simulator& ebos_simulator,
1971 std::vector<EvalWell>& mob,
1972 DeferredLogger& deferred_logger)
const
1974 const int cell_idx = this->well_cells_[perf];
1975 const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1976 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1980 if (this->isInjector()) {
1982 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1983 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1984 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
1987 if (PolymerModule::hasPlyshlog()) {
1990 if (this->isInjector() && this->wpolymer() == 0.) {
1995 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
1996 const EvalWell& bhp = this->getBhp();
1998 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
1999 double perf_dis_gas_rate = 0.;
2000 double perf_vap_oil_rate = 0.;
2001 double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
2002 const double Tw = this->well_index_[perf] * trans_mult;
2003 computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
2004 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
2006 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
2007 const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
2008 const auto& scaled_drainage_info =
2009 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
2010 const double swcr = scaled_drainage_info.Swcr;
2011 const EvalWell poro = this->extendEval(int_quant.porosity());
2012 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
2014 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
2015 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2016 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
2018 if (PolymerModule::hasShrate()) {
2021 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
2023 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
2024 int_quant.pvtRegionIndex(),
2027 mob[waterCompIdx] /= shear_factor;
2031 template<
typename TypeTag>
2033 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
2040 typename SparseMatrixAdapter::MatrixBlock tmpMat;
2041 Dune::DynamicMatrix<Scalar> tmp;
2042 for (
auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
2044 const auto row_index = colC.index();
2046 for (
auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
2048 Detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp);
2049 Detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
2050 jacobian.addToBlock( row_index, colB.index(), tmpMat );
2059 template<
typename TypeTag>
2060 typename StandardWell<TypeTag>::EvalWell
2061 StandardWell<TypeTag>::
2062 pskinwater(
const double throughput,
2063 const EvalWell& water_velocity,
2064 DeferredLogger& deferred_logger)
const
2066 if constexpr (Base::has_polymermw) {
2067 const int water_table_id = this->well_ecl_.getPolymerProperties().m_skprwattable;
2068 if (water_table_id <= 0) {
2069 OPM_DEFLOG_THROW(std::runtime_error,
"Unused SKPRWAT table id used for well " << name(), deferred_logger);
2071 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
2072 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2074 EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
2075 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
2078 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2079 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2087 template<
typename TypeTag>
2088 typename StandardWell<TypeTag>::EvalWell
2089 StandardWell<TypeTag>::
2090 pskin(
const double throughput,
2091 const EvalWell& water_velocity,
2092 const EvalWell& poly_inj_conc,
2093 DeferredLogger& deferred_logger)
const
2095 if constexpr (Base::has_polymermw) {
2096 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
2097 const EvalWell water_velocity_abs = abs(water_velocity);
2098 if (poly_inj_conc == 0.) {
2099 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
2101 const int polymer_table_id = this->well_ecl_.getPolymerProperties().m_skprpolytable;
2102 if (polymer_table_id <= 0) {
2103 OPM_DEFLOG_THROW(std::runtime_error,
"Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
2105 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
2106 const double reference_concentration = skprpolytable.refConcentration;
2107 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2109 EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
2110 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2111 if (poly_inj_conc == reference_concentration) {
2112 return sign * pskin_poly;
2115 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2116 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2117 return sign * pskin;
2119 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2120 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2128 template<
typename TypeTag>
2129 typename StandardWell<TypeTag>::EvalWell
2130 StandardWell<TypeTag>::
2131 wpolymermw(
const double throughput,
2132 const EvalWell& water_velocity,
2133 DeferredLogger& deferred_logger)
const
2135 if constexpr (Base::has_polymermw) {
2136 const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
2137 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2138 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2139 EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
2140 if (this->wpolymer() == 0.) {
2141 return molecular_weight;
2143 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2144 return molecular_weight;
2146 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2147 "while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
2155 template<
typename TypeTag>
2157 StandardWell<TypeTag>::
2158 updateWaterThroughput(
const double dt, WellState &well_state)
const
2160 if constexpr (Base::has_polymermw) {
2161 if (this->isInjector()) {
2162 auto& ws = well_state.well(this->index_of_well_);
2163 auto& perf_water_throughput = ws.perf_data.water_throughput;
2164 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2165 const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
2167 if (perf_water_vel > 0.) {
2168 perf_water_throughput[perf] += perf_water_vel * dt;
2179 template<
typename TypeTag>
2181 StandardWell<TypeTag>::
2182 handleInjectivityRate(
const Simulator& ebosSimulator,
2184 std::vector<EvalWell>& cq_s)
const
2186 const int cell_idx = this->well_cells_[perf];
2187 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2188 const auto& fs = int_quants.fluidState();
2189 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2190 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2191 const int wat_vel_index = Bhp + 1 + perf;
2192 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2196 cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
2202 template<
typename TypeTag>
2204 StandardWell<TypeTag>::
2205 handleInjectivityEquations(
const Simulator& ebosSimulator,
2206 const WellState& well_state,
2208 const EvalWell& water_flux_s,
2209 DeferredLogger& deferred_logger)
2211 const int cell_idx = this->well_cells_[perf];
2212 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2213 const auto& fs = int_quants.fluidState();
2214 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2215 const EvalWell water_flux_r = water_flux_s / b_w;
2216 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2217 const EvalWell water_velocity = water_flux_r / area;
2218 const int wat_vel_index = Bhp + 1 + perf;
2221 const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
2222 this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
2224 const auto& ws = well_state.well(this->index_of_well_);
2225 const auto& perf_data = ws.perf_data;
2226 const auto& perf_water_throughput = perf_data.water_throughput;
2227 const double throughput = perf_water_throughput[perf];
2228 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2230 EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
2231 poly_conc.setValue(this->wpolymer());
2234 const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
2235 - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
2237 this->resWell_[0][pskin_index] = eq_pskin.value();
2238 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
2239 this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
2240 this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
2244 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
2245 this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
2253 template<
typename TypeTag>
2255 StandardWell<TypeTag>::
2256 checkConvergenceExtraEqs(
const std::vector<double>& res,
2257 ConvergenceReport& report)
const
2262 if constexpr (Base::has_polymermw) {
2263 this->checkConvergencePolyMW(res, report, this->param_.max_residual_allowed_);
2271 template<
typename TypeTag>
2273 StandardWell<TypeTag>::
2274 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2275 const IntensiveQuantities& int_quants,
2276 const WellState& well_state,
2278 std::vector<RateVector>& connectionRates,
2279 DeferredLogger& deferred_logger)
const
2282 EvalWell cq_s_polymw = cq_s_poly;
2283 if (this->isInjector()) {
2284 const int wat_vel_index = Bhp + 1 + perf;
2285 const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
2286 if (water_velocity > 0.) {
2287 const auto& ws = well_state.well(this->index_of_well_);
2288 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2289 const double throughput = perf_water_throughput[perf];
2290 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2291 cq_s_polymw *= molecular_weight;
2297 }
else if (this->isProducer()) {
2298 if (cq_s_polymw < 0.) {
2299 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2306 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2314 template<
typename TypeTag>
2315 std::optional<double>
2316 StandardWell<TypeTag>::
2317 computeBhpAtThpLimitProd(
const WellState& well_state,
2318 const Simulator& ebos_simulator,
2319 const SummaryState& summary_state,
2320 DeferredLogger& deferred_logger)
const
2322 return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
2325 this->getALQ(well_state));
2328 template<
typename TypeTag>
2329 std::optional<double>
2330 StandardWell<TypeTag>::
2331 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
2332 const SummaryState& summary_state,
2333 DeferredLogger& deferred_logger,
2334 double alq_value)
const
2337 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2343 std::vector<double> rates(3);
2344 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2348 return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
2356 template<
typename TypeTag>
2357 std::optional<double>
2358 StandardWell<TypeTag>::
2359 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
2360 const SummaryState& summary_state,
2361 DeferredLogger& deferred_logger)
const
2364 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2370 std::vector<double> rates(3);
2371 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2375 return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitInj(frates,
2384 template<
typename TypeTag>
2386 StandardWell<TypeTag>::
2387 iterateWellEqWithControl(
const Simulator& ebosSimulator,
2389 const Well::InjectionControls& inj_controls,
2390 const Well::ProductionControls& prod_controls,
2391 WellState& well_state,
2392 const GroupState& group_state,
2393 DeferredLogger& deferred_logger)
2395 const int max_iter = this->param_.max_inner_iter_wells_;
2399 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2401 auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger);
2403 converged = report.converged();
2409 solveEqAndUpdateWellState(well_state, deferred_logger);
2416 initPrimaryVariablesEvaluation();
2417 }
while (it < max_iter);
2423 template<
typename TypeTag>
2430 std::vector<double> well_q_s(this->num_components_, 0.);
2431 const EvalWell& bhp = this->getBhp();
2432 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2433 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2434 const int cell_idx = this->well_cells_[perf];
2435 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2436 std::vector<Scalar> mob(this->num_components_, 0.);
2437 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
2438 std::vector<Scalar> cq_s(this->num_components_, 0.);
2439 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
2440 const double Tw = this->well_index_[perf] * trans_mult;
2441 computePerfRateScalar(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2442 cq_s, deferred_logger);
2443 for (
int comp = 0; comp < this->num_components_; ++comp) {
2444 well_q_s[comp] += cq_s[comp];
2447 const auto& comm = this->parallel_well_info_.communication();
2448 if (comm.size() > 1)
2450 comm.sum(well_q_s.data(), well_q_s.size());
2459 template <
typename TypeTag>
2463 const std::function<
double(
const double)>& connPICalc,
2464 const std::vector<EvalWell>& mobility,
2465 double* connPI)
const
2468 const int np = this->number_of_phases_;
2469 for (
int p = 0; p < np; ++p) {
2472 const auto connMob =
2473 mobility[ this->flowPhaseToEbosCompIdx(p) ].value()
2474 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
2476 connPI[p] = connPICalc(connMob);
2479 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2480 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2482 const auto io = pu.phase_pos[Oil];
2483 const auto ig = pu.phase_pos[Gas];
2485 const auto vapoil = connPI[ig] * fs.Rv().value();
2486 const auto disgas = connPI[io] * fs.Rs().value();
2488 connPI[io] += vapoil;
2489 connPI[ig] += disgas;
2497 template <
typename TypeTag>
2499 StandardWell<TypeTag>::
2500 computeConnLevelInjInd(
const typename StandardWell<TypeTag>::FluidState& fs,
2501 const Phase preferred_phase,
2502 const std::function<
double(
const double)>& connIICalc,
2503 const std::vector<EvalWell>& mobility,
2505 DeferredLogger& deferred_logger)
const
2511 if (preferred_phase == Phase::GAS) {
2512 phase_pos = pu.phase_pos[Gas];
2514 else if (preferred_phase == Phase::OIL) {
2515 phase_pos = pu.phase_pos[Oil];
2517 else if (preferred_phase == Phase::WATER) {
2518 phase_pos = pu.phase_pos[Water];
2521 OPM_DEFLOG_THROW(NotImplemented,
2522 "Unsupported Injector Type ("
2523 <<
static_cast<int>(preferred_phase)
2524 <<
") for well " << this->name()
2525 <<
" during connection I.I. calculation",
2529 const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
2530 const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
2531 connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition: ConvergenceReport.hpp:36
Definition: DeferredLogger.hpp:57
Definition: StandardWell.hpp:65
Collect per-connection static information to enable calculating connection-level or well-level produc...
Definition: WellProdIndexCalculator.hpp:36
double connectionProdIndStandard(const std::size_t connIdx, const double connMobility) const
Compute connection-level steady-state productivity index value using dynamic phase mobility.
Definition: WellProdIndexCalculator.cpp:106
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33