22#include <opm/common/utility/numeric/RootFinders.hpp>
23#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
25#include <opm/simulators/wells/VFPHelpers.hpp>
34 template<
typename TypeTag>
35 StandardWell<TypeTag>::
36 StandardWell(
const Well& well,
37 const ParallelWellInfo& pw_info,
39 const ModelParameters& param,
40 const RateConverterType& rate_converter,
41 const int pvtRegionIdx,
42 const int num_components,
44 const int index_of_well,
45 const std::vector<PerforationData>& perf_data)
46 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
47 , 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,
65 const bool changed_to_open_this_step)
67 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
68 this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
75 template<
typename TypeTag>
76 void StandardWell<TypeTag>::
77 initPrimaryVariablesEvaluation()
const
79 this->StdWellEval::initPrimaryVariablesEvaluation();
86 template<
typename TypeTag>
88 StandardWell<TypeTag>::
89 computePerfRateEval(
const IntensiveQuantities& intQuants,
90 const std::vector<EvalWell>& mob,
95 std::vector<EvalWell>& cq_s,
96 double& perf_dis_gas_rate,
97 double& perf_vap_oil_rate,
98 double& perf_vap_wat_rate,
99 DeferredLogger& deferred_logger)
const
101 const auto& fs = intQuants.fluidState();
102 const EvalWell pressure = this->extendEval(this->getPerfCellPressure(fs));
103 const EvalWell rs = this->extendEval(fs.Rs());
104 const EvalWell rv = this->extendEval(fs.Rv());
105 const EvalWell rvw = this->extendEval(fs.Rvw());
107 std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
108 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
109 if (!FluidSystem::phaseIsActive(phaseIdx)) {
112 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
113 b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
115 if constexpr (has_solvent) {
116 b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
119 if constexpr (has_zFraction) {
120 if (this->isInjector()) {
121 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
122 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
123 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
127 EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
129 if (this->isInjector()) {
130 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
131 skin_pressure = this->primary_variables_evaluation_[pskin_index];
136 std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
137 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
138 cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
160 template<
typename TypeTag>
162 StandardWell<TypeTag>::
163 computePerfRateScalar(
const IntensiveQuantities& intQuants,
164 const std::vector<Scalar>& mob,
169 std::vector<Scalar>& cq_s,
170 DeferredLogger& deferred_logger)
const
172 const auto& fs = intQuants.fluidState();
173 const Scalar pressure = this->getPerfCellPressure(fs).value();
174 const Scalar rs = fs.Rs().value();
175 const Scalar rv = fs.Rv().value();
176 const Scalar rvw = fs.Rvw().value();
177 std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
178 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
179 if (!FluidSystem::phaseIsActive(phaseIdx)) {
182 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
183 b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
185 if constexpr (has_solvent) {
186 b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
189 if constexpr (has_zFraction) {
190 if (this->isInjector()) {
191 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
192 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
193 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
197 Scalar skin_pressure =0.0;
199 if (this->isInjector()) {
200 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
201 skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
205 Scalar perf_dis_gas_rate = 0.0;
206 Scalar perf_vap_oil_rate = 0.0;
207 Scalar perf_vap_wat_rate = 0.0;
210 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
211 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
212 cmix_s[componentIdx] = getValue(this->wellSurfaceVolumeFraction(componentIdx));
234 template<
typename TypeTag>
235 template<
class Value>
237 StandardWell<TypeTag>::
238 computePerfRate(
const std::vector<Value>& mob,
239 const Value& pressure,
244 std::vector<Value>& b_perfcells_dense,
248 const Value& skin_pressure,
249 const std::vector<Value>& cmix_s,
250 std::vector<Value>& cq_s,
251 double& perf_dis_gas_rate,
252 double& perf_vap_oil_rate,
253 double& perf_vap_wat_rate,
254 DeferredLogger& deferred_logger)
const
257 const Value well_pressure = bhp + this->perf_pressure_diffs_[perf];
258 Value drawdown = pressure - well_pressure;
259 if (this->isInjector()) {
260 drawdown += skin_pressure;
264 if ( drawdown > 0 ) {
266 if (!allow_cf && this->isInjector()) {
271 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
272 const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
273 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
276 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
277 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
278 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
279 const Value cq_sOil = cq_s[oilCompIdx];
280 const Value cq_sGas = cq_s[gasCompIdx];
281 const Value dis_gas = rs * cq_sOil;
282 const Value vap_oil = rv * cq_sGas;
284 cq_s[gasCompIdx] += dis_gas;
285 cq_s[oilCompIdx] += vap_oil;
288 if (this->isProducer()) {
289 perf_dis_gas_rate = getValue(dis_gas);
290 perf_vap_oil_rate = getValue(vap_oil);
292 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
293 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
294 const Value vap_wat = rvw * cq_sGas;
295 cq_s[waterCompIdx] += vap_wat;
296 if (this->isProducer())
297 perf_vap_wat_rate = getValue(vap_wat);
303 if (!allow_cf && this->isProducer()) {
308 Value total_mob_dense = mob[0];
309 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
310 total_mob_dense += mob[componentIdx];
314 const Value cqt_i = - Tw * (total_mob_dense * drawdown);
317 Value volumeRatio = bhp * 0.0;
319 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
320 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
321 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
324 if constexpr (Indices::enableSolvent) {
325 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
328 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
329 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
330 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
332 const Value d = 1.0 - rv * rs;
335 std::ostringstream sstr;
336 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
337 <<
" during computePerfRate calculations with rs " << rs
338 <<
", rv " << rv <<
" and pressure " << pressure
339 <<
" obtaining d " << d
340 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
341 <<
" for this connection.";
342 deferred_logger.debug(sstr.str());
344 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
345 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
347 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
348 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
351 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
352 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
353 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
355 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
356 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
357 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
362 Value cqt_is = cqt_i/volumeRatio;
363 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
364 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
368 if (this->isProducer()) {
369 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
370 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
371 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
380 const double d = 1.0 - getValue(rv) * getValue(rs);
383 std::ostringstream sstr;
384 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
385 <<
" during computePerfRate calculations with rs " << rs
386 <<
", rv " << rv <<
" and pressure " << pressure
387 <<
" obtaining d " << d
388 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
389 <<
" for this connection.";
390 deferred_logger.debug(sstr.str());
394 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
397 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
399 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
404 perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
407 else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
409 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
410 perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]);
417 template<
typename TypeTag>
419 StandardWell<TypeTag>::
420 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
422 const Well::InjectionControls& ,
423 const Well::ProductionControls& ,
424 WellState& well_state,
425 const GroupState& group_state,
426 DeferredLogger& deferred_logger)
430 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
436 this->resWell_ = 0.0;
438 assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
444 template<
typename TypeTag>
446 StandardWell<TypeTag>::
447 assembleWellEqWithoutIterationImpl(
const Simulator& ebosSimulator,
449 WellState& well_state,
450 const GroupState& group_state,
451 DeferredLogger& deferred_logger)
454 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
455 const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
457 auto& ws = well_state.well(this->index_of_well_);
459 ws.vaporized_oil_rate = 0;
460 ws.dissolved_gas_rate = 0;
461 ws.vaporized_wat_rate = 0;
463 const int np = this->number_of_phases_;
465 std::vector<RateVector> connectionRates = this->connectionRates_;
466 auto& perf_data = ws.perf_data;
467 auto& perf_rates = perf_data.phase_rates;
468 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
470 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
471 EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
472 EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
473 calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
476 if constexpr (has_polymer && Base::has_polymermw) {
477 if (this->isInjector()) {
478 handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
481 const int cell_idx = this->well_cells_[perf];
482 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
484 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
486 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
489 this->resWell_[0][componentIdx] += cq_s_effective.value();
492 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
494 this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq);
495 this->duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
498 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
499 this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
503 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
504 auto& perf_rate_solvent = perf_data.solvent_rates;
505 perf_rate_solvent[perf] = cq_s[componentIdx].value();
507 perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
511 if constexpr (has_zFraction) {
512 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
513 this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
518 this->connectionRates_ = connectionRates;
523 const auto& comm = this->parallel_well_info_.communication();
524 ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
525 ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
526 ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
530 wellhelpers::sumDistributedWellEntries(this->duneD_[0][0], this->resWell_[0],
531 this->parallel_well_info_.communication());
533 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
536 EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
537 if (FluidSystem::numActivePhases() > 1) {
539 resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
541 resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
542 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
543 this->duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
545 this->resWell_[0][componentIdx] += resWell_loc.value();
548 const auto& summaryState = ebosSimulator.vanguard().summaryState();
549 const Schedule& schedule = ebosSimulator.vanguard().schedule();
550 this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger);
555 this->invDuneD_ = this->duneD_;
556 detail::invertMatrix(this->invDuneD_[0][0]);
557 }
catch (NumericalProblem&) {
559 this->invDuneD_[0][0] = 0.0;
560 for (
size_t i = 0; i < this->invDuneD_[0][0].rows(); ++i) {
561 this->invDuneD_[0][0][i][i] = 1.0;
564 OPM_DEFLOG_THROW(NumericalIssue,
"Error when inverting local well equations for well " + name(), deferred_logger);
571 template<
typename TypeTag>
573 StandardWell<TypeTag>::
574 calculateSinglePerf(
const Simulator& ebosSimulator,
576 WellState& well_state,
577 std::vector<RateVector>& connectionRates,
578 std::vector<EvalWell>& cq_s,
579 EvalWell& water_flux_s,
580 EvalWell& cq_s_zfrac_effective,
581 DeferredLogger& deferred_logger)
const
583 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
584 const EvalWell& bhp = this->getBhp();
585 const int cell_idx = this->well_cells_[perf];
586 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
587 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
588 getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
590 double perf_dis_gas_rate = 0.;
591 double perf_vap_oil_rate = 0.;
592 double perf_vap_wat_rate = 0.;
593 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
594 const double Tw = this->well_index_[perf] * trans_mult;
595 computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
596 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
598 auto& ws = well_state.well(this->index_of_well_);
599 auto& perf_data = ws.perf_data;
600 if constexpr (has_polymer && Base::has_polymermw) {
601 if (this->isInjector()) {
604 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
605 water_flux_s = cq_s[water_comp_idx];
608 handleInjectivityRate(ebosSimulator, perf, cq_s);
613 if (this->isProducer()) {
614 ws.dissolved_gas_rate += perf_dis_gas_rate;
615 ws.vaporized_oil_rate += perf_vap_oil_rate;
616 ws.vaporized_wat_rate += perf_vap_wat_rate;
619 if constexpr (has_energy) {
620 connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
623 if constexpr (has_energy) {
625 auto fs = intQuants.fluidState();
626 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
627 if (!FluidSystem::phaseIsActive(phaseIdx)) {
632 EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
633 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
634 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
635 if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
636 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
639 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
640 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
645 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
647 std::ostringstream sstr;
648 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
649 <<
" during calculateSinglePerf with rs " << fs.Rs()
650 <<
", rv " << fs.Rv()
651 <<
" obtaining d " << d
652 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
653 <<
" for this connection.";
654 deferred_logger.debug(sstr.str());
655 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
657 if(FluidSystem::gasPhaseIdx == phaseIdx) {
658 cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
659 }
else if(FluidSystem::oilPhaseIdx == phaseIdx) {
661 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
667 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
669 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
670 fs.setTemperature(this->well_ecl_.temperature());
671 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
672 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
673 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
674 paramCache.setRegionIndex(pvtRegionIdx);
675 paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
676 paramCache.updatePhase(fs, phaseIdx);
678 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
679 fs.setDensity(phaseIdx, rho);
680 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
681 fs.setEnthalpy(phaseIdx, h);
684 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
685 connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
689 if constexpr (has_polymer) {
691 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
692 EvalWell cq_s_poly = cq_s[waterCompIdx];
693 if (this->isInjector()) {
694 cq_s_poly *= this->wpolymer();
696 cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
699 auto& perf_rate_polymer = perf_data.polymer_rates;
700 perf_rate_polymer[perf] = cq_s_poly.value();
702 cq_s_poly *= this->well_efficiency_factor_;
703 connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
705 if constexpr (Base::has_polymermw) {
706 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
710 if constexpr (has_foam) {
712 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
713 EvalWell cq_s_foam = cq_s[gasCompIdx] * this->well_efficiency_factor_;
714 if (this->isInjector()) {
715 cq_s_foam *= this->wfoam();
717 cq_s_foam *= this->extendEval(intQuants.foamConcentration());
719 connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
722 if constexpr (has_zFraction) {
724 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
725 cq_s_zfrac_effective = cq_s[gasCompIdx];
726 if (this->isInjector()) {
727 cq_s_zfrac_effective *= this->wsolvent();
728 }
else if (cq_s_zfrac_effective.value() != 0.0) {
729 const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
730 cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
732 auto& perf_rate_solvent = perf_data.solvent_rates;
733 perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
735 cq_s_zfrac_effective *= this->well_efficiency_factor_;
736 connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
739 if constexpr (has_brine) {
741 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
743 EvalWell cq_s_sm = cq_s[waterCompIdx] - perf_vap_wat_rate;
744 if (this->isInjector()) {
745 cq_s_sm *= this->wsalt();
747 cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration());
750 auto& perf_rate_brine = perf_data.brine_rates;
751 perf_rate_brine[perf] = cq_s_sm.value();
753 cq_s_sm *= this->well_efficiency_factor_;
754 connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
757 if constexpr (has_micp) {
758 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
759 EvalWell cq_s_microbe = cq_s[waterCompIdx];
760 if (this->isInjector()) {
761 cq_s_microbe *= this->wmicrobes();
763 cq_s_microbe *= this->extendEval(intQuants.microbialConcentration());
765 connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe);
766 EvalWell cq_s_oxygen = cq_s[waterCompIdx];
767 if (this->isInjector()) {
768 cq_s_oxygen *= this->woxygen();
770 cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration());
772 connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen);
773 EvalWell cq_s_urea = cq_s[waterCompIdx];
774 if (this->isInjector()) {
775 cq_s_urea *= this->wurea();
777 cq_s_urea *= this->extendEval(intQuants.ureaConcentration());
779 connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea);
783 perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf];
789 template<
typename TypeTag>
791 StandardWell<TypeTag>::
792 getMobilityEval(
const Simulator& ebosSimulator,
794 std::vector<EvalWell>& mob,
795 DeferredLogger& deferred_logger)
const
797 const int cell_idx = this->well_cells_[perf];
798 assert (
int(mob.size()) == this->num_components_);
799 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
800 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
804 const int satid = this->saturation_table_number_[perf] - 1;
805 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
806 if( satid == satid_elem ) {
808 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
809 if (!FluidSystem::phaseIsActive(phaseIdx)) {
813 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
814 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
817 mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
821 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
822 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
823 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
826 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
829 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
830 if (!FluidSystem::phaseIsActive(phaseIdx)) {
834 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
835 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
839 if constexpr (has_solvent) {
840 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
845 if constexpr (has_polymer) {
846 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
847 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
852 if constexpr (!Base::has_polymermw) {
853 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
858 template<
typename TypeTag>
860 StandardWell<TypeTag>::
861 getMobilityScalar(
const Simulator& ebosSimulator,
863 std::vector<Scalar>& mob,
864 DeferredLogger& deferred_logger)
const
866 const int cell_idx = this->well_cells_[perf];
867 assert (
int(mob.size()) == this->num_components_);
868 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
869 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
873 const int satid = this->saturation_table_number_[perf] - 1;
874 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
875 if( satid == satid_elem ) {
877 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
878 if (!FluidSystem::phaseIsActive(phaseIdx)) {
882 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
883 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
886 mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
890 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
891 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
892 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
895 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
898 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
899 if (!FluidSystem::phaseIsActive(phaseIdx)) {
903 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
904 mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx));
908 if constexpr (has_solvent) {
909 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
914 if constexpr (has_polymer) {
915 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
916 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
921 if constexpr (!Base::has_polymermw) {
922 std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
923 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
924 for (
size_t i = 0; i < mob.size(); ++i) {
925 mob[i] = getValue(mob_eval[i]);
933 template<
typename TypeTag>
935 StandardWell<TypeTag>::
936 updateWellState(
const BVectorWell& dwells,
937 WellState& well_state,
938 DeferredLogger& deferred_logger)
const
940 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
942 updatePrimaryVariablesNewton(dwells, well_state, deferred_logger);
944 updateWellStateFromPrimaryVariables(well_state, deferred_logger);
945 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
952 template<
typename TypeTag>
954 StandardWell<TypeTag>::
955 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
957 DeferredLogger& deferred_logger)
const
959 const double dFLimit = this->param_.dwell_fraction_max_;
960 const double dBHPLimit = this->param_.dbhp_max_rel_;
961 this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
963 updateExtraPrimaryVariables(dwells);
965 for (
double v : this->primary_variables_) {
967 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after newton update well: " << this->name(), deferred_logger);
976 template<
typename TypeTag>
978 StandardWell<TypeTag>::
979 updateExtraPrimaryVariables(
const BVectorWell& dwells)
const
982 if constexpr (Base::has_polymermw) {
983 this->updatePrimaryVariablesPolyMW(dwells);
991 template<
typename TypeTag>
993 StandardWell<TypeTag>::
994 updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger)
const
996 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger);
999 if constexpr (Base::has_polymermw) {
1000 this->updateWellStateFromPrimaryVariablesPolyMW(well_state);
1008 template<
typename TypeTag>
1010 StandardWell<TypeTag>::
1011 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1016 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1017 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1019 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1020 std::vector<Scalar> mob(this->num_components_, 0.0);
1021 getMobilityScalar(ebos_simulator, perf, mob, deferred_logger);
1023 const int cell_idx = this->well_cells_[perf];
1024 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1025 const auto& fs = int_quantities.fluidState();
1027 double p_r = this->getPerfCellPressure(fs).value();
1030 std::vector<double> b_perf(this->num_components_);
1031 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1032 if (!FluidSystem::phaseIsActive(phase)) {
1035 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1036 b_perf[comp_idx] = fs.invB(phase).value();
1038 if constexpr (has_solvent) {
1039 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
1043 const double h_perf = this->perf_pressure_diffs_[perf];
1044 const double pressure_diff = p_r - h_perf;
1049 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1050 deferred_logger.debug(
"CROSSFLOW_IPR",
1051 "cross flow found when updateIPR for well " + name()
1052 +
" . The connection is ignored in IPR calculations");
1058 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1060 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1061 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1062 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1063 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1064 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1065 ipr_b_perf[comp_idx] += tw_mob;
1069 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1070 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1071 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1072 const double rs = (fs.Rs()).value();
1073 const double rv = (fs.Rv()).value();
1075 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1076 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1078 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1079 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1081 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1082 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1084 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1085 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1088 for (
size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1089 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1090 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1093 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1094 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1098 template<
typename TypeTag>
1100 StandardWell<TypeTag>::
1101 checkOperabilityUnderBHPLimit(
const WellState& well_state,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1103 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1104 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1107 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1108 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1111 double total_ipr_mass_rate = 0.0;
1112 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1114 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1118 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1119 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1121 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1122 total_ipr_mass_rate += ipr_rate * rho;
1124 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1125 this->operability_status_.operable_under_only_bhp_limit =
false;
1129 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1133 std::vector<double> well_rates_bhp_limit;
1134 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1136 this->adaptRatesForVFP(well_rates_bhp_limit);
1137 const double thp = this->calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
1138 const double thp_limit = this->getTHPConstraint(summaryState);
1139 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1140 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1151 this->operability_status_.operable_under_only_bhp_limit =
true;
1152 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1160 template<
typename TypeTag>
1162 StandardWell<TypeTag>::
1163 checkOperabilityUnderTHPLimit(
const Simulator& ebos_simulator,
const WellState& well_state, DeferredLogger& deferred_logger)
1165 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1166 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger)
1167 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1170 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1172 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1173 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1175 const double thp_limit = this->getTHPConstraint(summaryState);
1176 if (this->isProducer() && *obtain_bhp < thp_limit) {
1177 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1178 +
" bars is SMALLER than thp limit "
1179 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1180 +
" bars as a producer for well " + name();
1181 deferred_logger.debug(msg);
1183 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1184 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1185 +
" bars is LARGER than thp limit "
1186 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1187 +
" bars as a injector for well " + name();
1188 deferred_logger.debug(msg);
1191 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1192 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1193 if (!this->wellIsStopped()) {
1194 const double thp_limit = this->getTHPConstraint(summaryState);
1195 deferred_logger.debug(
" could not find bhp value at thp limit "
1196 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1197 +
" bar for well " + name() +
", the well might need to be closed ");
1206 template<
typename TypeTag>
1208 StandardWell<TypeTag>::
1209 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1211 bool all_drawdown_wrong_direction =
true;
1213 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1214 const int cell_idx = this->well_cells_[perf];
1215 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1216 const auto& fs = intQuants.fluidState();
1218 const double pressure = this->getPerfCellPressure(fs).value();
1219 const double bhp = this->getBhp().value();
1222 const double well_pressure = bhp + this->perf_pressure_diffs_[perf];
1223 const double drawdown = pressure - well_pressure;
1228 if ( (drawdown < 0. && this->isInjector()) ||
1229 (drawdown > 0. && this->isProducer()) ) {
1230 all_drawdown_wrong_direction =
false;
1235 const auto& comm = this->parallel_well_info_.communication();
1236 if (comm.size() > 1)
1238 all_drawdown_wrong_direction =
1239 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1242 return all_drawdown_wrong_direction;
1248 template<
typename TypeTag>
1250 StandardWell<TypeTag>::
1251 canProduceInjectWithCurrentBhp(
const Simulator& ebos_simulator,
1252 const WellState& well_state,
1253 DeferredLogger& deferred_logger)
1255 const double bhp = well_state.well(this->index_of_well_).bhp;
1256 std::vector<double> well_rates;
1257 computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
1259 const double sign = (this->isProducer()) ? -1. : 1.;
1260 const double threshold = sign * std::numeric_limits<double>::min();
1262 bool can_produce_inject =
false;
1263 for (
const auto value : well_rates) {
1264 if (this->isProducer() && value < threshold) {
1265 can_produce_inject =
true;
1267 }
else if (this->isInjector() && value > threshold) {
1268 can_produce_inject =
true;
1273 if (!can_produce_inject) {
1274 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1277 return can_produce_inject;
1284 template<
typename TypeTag>
1286 StandardWell<TypeTag>::
1287 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1289 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1295 template<
typename TypeTag>
1297 StandardWell<TypeTag>::
1298 computePropertiesForWellConnectionPressures(
const Simulator& ebosSimulator,
1299 const WellState& well_state,
1300 std::vector<double>& b_perf,
1301 std::vector<double>& rsmax_perf,
1302 std::vector<double>& rvmax_perf,
1303 std::vector<double>& rvwmax_perf,
1304 std::vector<double>& surf_dens_perf)
const
1306 const int nperf = this->number_of_perforations_;
1308 b_perf.resize(nperf * this->num_components_);
1309 surf_dens_perf.resize(nperf * this->num_components_);
1310 const auto& ws = well_state.well(this->index_of_well_);
1312 const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1313 const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1314 const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1317 if (oilPresent && gasPresent) {
1318 rsmax_perf.resize(nperf);
1319 rvmax_perf.resize(nperf);
1322 if (waterPresent && gasPresent) {
1323 rvwmax_perf.resize(nperf);
1327 const auto& perf_press = ws.perf_data.pressure;
1328 auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp,
1332 for (
int perf = 0; perf < nperf; ++perf) {
1333 const int cell_idx = this->well_cells_[perf];
1334 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1335 const auto& fs = intQuants.fluidState();
1337 const double p_avg = (perf_press[perf] + p_above[perf])/2;
1338 const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1339 const double saltConcentration = fs.saltConcentration().value();
1342 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1343 b_perf[ waterCompIdx + perf * this->num_components_] =
1344 FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration);
1348 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1349 const int gaspos = gasCompIdx + perf * this->num_components_;
1351 if (oilPresent && waterPresent) {
1352 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1353 const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]);
1354 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1355 rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1359 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1361 rv = oilrate / gasrate;
1363 rv = std::min(rv, rvmax_perf[perf]);
1365 if (waterrate > 0) {
1366 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1368 rvw = waterrate / gasrate;
1370 rvw = std::min(rvw, rvwmax_perf[perf]);
1372 if (rv > 0.0 || rvw > 0.0){
1373 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, rvw);
1376 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1378 }
else if (oilPresent) {
1380 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1381 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1383 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1386 rv = oilrate / gasrate;
1388 rv = std::min(rv, rvmax_perf[perf]);
1390 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, 0.0 );
1393 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1395 }
else if (waterPresent) {
1397 const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]);
1398 rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1399 if (waterrate > 0) {
1400 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1403 rvw = waterrate / gasrate;
1405 rvw = std::min(rvw, rvwmax_perf[perf]);
1407 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, 0.0 , rvw);
1410 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1414 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1419 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1420 const int oilpos = oilCompIdx + perf * this->num_components_;
1422 rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1423 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1425 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1428 rs = gasrate / oilrate;
1430 rs = std::min(rs, rsmax_perf[perf]);
1431 b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1433 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1436 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1441 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1442 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1446 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1447 surf_dens_perf[this->num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() );
1451 if constexpr (has_solvent) {
1452 b_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1453 surf_dens_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity();
1462 template<
typename TypeTag>
1466 const std::vector<double>& B_avg,
1468 const bool relax_tolerance)
const
1472 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1474 std::vector<double> res;
1477 this->param_.max_residual_allowed_,
1478 this->param_.tolerance_wells_,
1479 this->param_.relaxed_tolerance_flow_well_,
1483 checkConvergenceExtraEqs(res, report);
1492 template<
typename TypeTag>
1500 auto fluidState = [&ebosSimulator,
this](
const int perf)
1502 const auto cell_idx = this->well_cells_[perf];
1503 return ebosSimulator.model()
1504 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
1507 const int np = this->number_of_phases_;
1508 auto setToZero = [np](
double* x) ->
void
1510 std::fill_n(x, np, 0.0);
1513 auto addVector = [np](
const double* src,
double* dest) ->
void
1515 std::transform(src, src + np, dest, dest, std::plus<>{});
1518 auto& ws = well_state.well(this->index_of_well_);
1519 auto& perf_data = ws.perf_data;
1520 auto* wellPI = ws.productivity_index.data();
1521 auto* connPI = perf_data.prod_index.data();
1525 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1526 auto subsetPerfID = 0;
1528 for (
const auto& perf : *this->perf_data_) {
1529 auto allPerfID = perf.ecl_index;
1531 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1536 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
1537 getMobilityEval(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1539 const auto& fs = fluidState(subsetPerfID);
1542 if (this->isInjector()) {
1543 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1544 mob, connPI, deferred_logger);
1547 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1550 addVector(connPI, wellPI);
1557 const auto& comm = this->parallel_well_info_.communication();
1558 if (comm.size() > 1) {
1559 comm.sum(wellPI, np);
1562 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1563 "Internal logic error in processing connections for PI/II");
1568 template<
typename TypeTag>
1570 StandardWell<TypeTag>::
1571 computeWellConnectionDensitesPressures(
const Simulator& ebosSimulator,
1572 const WellState& well_state,
1573 const std::vector<double>& b_perf,
1574 const std::vector<double>& rsmax_perf,
1575 const std::vector<double>& rvmax_perf,
1576 const std::vector<double>& rvwmax_perf,
1577 const std::vector<double>& surf_dens_perf,
1578 DeferredLogger& deferred_logger)
1581 const int nperf = this->number_of_perforations_;
1582 const int np = this->number_of_phases_;
1583 std::vector<double> perfRates(b_perf.size(),0.0);
1584 const auto& ws = well_state.well(this->index_of_well_);
1585 const auto& perf_data = ws.perf_data;
1586 const auto& perf_rates_state = perf_data.phase_rates;
1588 for (
int perf = 0; perf < nperf; ++perf) {
1589 for (
int comp = 0; comp < np; ++comp) {
1590 perfRates[perf * this->num_components_ + comp] = perf_rates_state[perf * np + this->ebosCompIdxToFlowCompIdx(comp)];
1594 if constexpr (has_solvent) {
1595 const auto& solvent_perf_rates_state = perf_data.solvent_rates;
1596 for (
int perf = 0; perf < nperf; ++perf) {
1597 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = solvent_perf_rates_state[perf];
1604 bool all_zero = std::all_of(perfRates.begin(), perfRates.end(), [](
double val) { return val == 0.0; });
1605 const auto& comm = this->parallel_well_info_.communication();
1606 if (comm.size() > 1)
1608 all_zero = (comm.min(all_zero ? 1 : 0) == 1);
1611 if ( all_zero && this->isProducer() ) {
1612 double total_tw = 0;
1613 for (
int perf = 0; perf < nperf; ++perf) {
1614 total_tw += this->well_index_[perf];
1616 if (comm.size() > 1)
1618 total_tw = comm.sum(total_tw);
1620 for (
int perf = 0; perf < nperf; ++perf) {
1621 const int cell_idx = this->well_cells_[perf];
1622 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1623 const auto& fs = intQuants.fluidState();
1624 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1625 double total_mobility = 0.0;
1626 for (
int p = 0; p < np; ++p) {
1627 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1628 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1630 if constexpr (has_solvent) {
1631 total_mobility += intQuants.solventInverseFormationVolumeFactor().value() * intQuants.solventMobility().value();
1633 for (
int p = 0; p < np; ++p) {
1634 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1635 perfRates[perf * this->num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1637 if constexpr (has_solvent) {
1638 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility;
1643 this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
1645 this->computeConnectionPressureDelta();
1652 template<
typename TypeTag>
1654 StandardWell<TypeTag>::
1655 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1656 const WellState& well_state,
1657 DeferredLogger& deferred_logger)
1662 std::vector<double> b_perf;
1663 std::vector<double> rsmax_perf;
1664 std::vector<double> rvmax_perf;
1665 std::vector<double> rvwmax_perf;
1666 std::vector<double> surf_dens_perf;
1667 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf);
1668 computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
1675 template<
typename TypeTag>
1677 StandardWell<TypeTag>::
1678 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
1680 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1684 BVectorWell dx_well(1);
1685 dx_well[0].resize(this->numWellEq_);
1686 this->invDuneD_.mv(this->resWell_, dx_well);
1688 updateWellState(dx_well, well_state, deferred_logger);
1695 template<
typename TypeTag>
1697 StandardWell<TypeTag>::
1698 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1699 const WellState& well_state,
1700 DeferredLogger& deferred_logger)
1702 updatePrimaryVariables(well_state, deferred_logger);
1703 initPrimaryVariablesEvaluation();
1704 computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
1705 this->computeAccumWell();
1710 template<
typename TypeTag>
1713 apply(
const BVector& x, BVector& Ax)
const
1715 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1717 if (this->param_.matrix_add_well_contributions_)
1722 assert( this->Bx_.size() == this->duneB_.N() );
1723 assert( this->invDrw_.size() == this->invDuneD_.N() );
1726 this->parallelB_.mv(x, this->Bx_);
1731 BVectorWell& invDBx = this->invDrw_;
1732 this->invDuneD_.mv(this->Bx_, invDBx);
1735 this->duneC_.mmtv(invDBx,Ax);
1741 template<
typename TypeTag>
1744 apply(BVector& r)
const
1746 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1748 assert( this->invDrw_.size() == this->invDuneD_.N() );
1751 this->invDuneD_.mv(this->resWell_, this->invDrw_);
1753 this->duneC_.mmtv(this->invDrw_, r);
1756 template<
typename TypeTag>
1761 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1763 BVectorWell resWell = this->resWell_;
1765 this->parallelB_.mmv(x, resWell);
1767 this->invDuneD_.mv(resWell, xw);
1774 template<
typename TypeTag>
1781 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1784 xw[0].resize(this->numWellEq_);
1786 recoverSolutionWell(x, xw);
1787 updateWellState(xw, well_state, deferred_logger);
1793 template<
typename TypeTag>
1798 std::vector<double>& well_flux,
1802 const int np = this->number_of_phases_;
1803 well_flux.resize(np, 0.0);
1805 const bool allow_cf = this->getAllowCrossFlow();
1807 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1808 const int cell_idx = this->well_cells_[perf];
1809 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1811 std::vector<Scalar> mob(this->num_components_, 0.);
1812 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
1813 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
1814 const double Tw = this->well_index_[perf] * trans_mult;
1816 std::vector<Scalar> cq_s(this->num_components_, 0.);
1817 computePerfRateScalar(intQuants, mob, bhp, Tw, perf, allow_cf,
1818 cq_s, deferred_logger);
1820 for(
int p = 0; p < np; ++p) {
1821 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
1824 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1829 template<
typename TypeTag>
1831 StandardWell<TypeTag>::
1832 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
1834 std::vector<double>& well_flux,
1835 DeferredLogger& deferred_logger)
const
1841 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1842 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1843 auto& ws = well_state_copy.well(this->index_of_well_);
1846 if (this->well_ecl_.isInjector()) {
1847 ws.injection_cmode = Well::InjectorCMode::BHP;
1849 ws.production_cmode = Well::ProducerCMode::BHP;
1854 const int np = this->number_of_phases_;
1855 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1856 for (
int phase = 0; phase < np; ++phase){
1857 well_state_copy.wellRates(this->index_of_well_)[phase]
1858 = sign * ws.well_potentials[phase];
1862 StandardWell<TypeTag> well(*
this);
1863 well.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1865 const double dt = ebosSimulator.timeStepSize();
1866 bool converged = well.iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
1868 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1869 " potentials are computed based on unconverged solution";
1870 deferred_logger.debug(msg);
1872 well.updatePrimaryVariables(well_state_copy, deferred_logger);
1873 well.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
1874 well.initPrimaryVariablesEvaluation();
1875 well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1881 template<
typename TypeTag>
1883 StandardWell<TypeTag>::
1884 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
1885 DeferredLogger& deferred_logger,
1886 const WellState &well_state)
const
1888 std::vector<double> potentials(this->number_of_phases_, 0.0);
1889 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1891 const auto& well = this->well_ecl_;
1892 if (well.isInjector()){
1893 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1894 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1895 if (bhp_at_thp_limit) {
1896 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1897 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1899 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1900 "Failed in getting converged thp based potential calculation for well "
1901 + name() +
". Instead the bhp based value is used");
1902 const double bhp = controls.bhp_limit;
1903 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1906 computeWellRatesWithThpAlqProd(
1907 ebos_simulator, summary_state,
1908 deferred_logger, potentials, this->getALQ(well_state)
1915 template<
typename TypeTag>
1917 StandardWell<TypeTag>::
1918 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &ebos_simulator,
1919 const SummaryState &summary_state,
1920 DeferredLogger &deferred_logger,
1921 std::vector<double> &potentials,
1925 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1926 ebos_simulator, summary_state, deferred_logger, alq);
1927 if (bhp_at_thp_limit) {
1928 const auto& controls = this->well_ecl_.productionControls(summary_state);
1929 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1930 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1933 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1934 "Failed in getting converged thp based potential calculation for well "
1935 + name() +
". Instead the bhp based value is used");
1936 const auto& controls = this->well_ecl_.productionControls(summary_state);
1937 bhp = controls.bhp_limit;
1938 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1943 template<
typename TypeTag>
1945 StandardWell<TypeTag>::
1946 computeWellRatesWithThpAlqProd(
const Simulator &ebos_simulator,
1947 const SummaryState &summary_state,
1948 DeferredLogger &deferred_logger,
1949 std::vector<double> &potentials,
1953 computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1960 template<
typename TypeTag>
1965 std::vector<double>& well_potentials,
1968 const int np = this->number_of_phases_;
1969 well_potentials.resize(np, 0.0);
1971 if (this->wellIsStopped()) {
1975 this->operability_status_.has_negative_potentials =
false;
1977 bool thp_controlled_well =
false;
1978 bool bhp_controlled_well =
false;
1979 const auto& ws = well_state.well(this->index_of_well_);
1980 if (this->isInjector()) {
1981 const Well::InjectorCMode& current = ws.injection_cmode;
1982 if (current == Well::InjectorCMode::THP) {
1983 thp_controlled_well =
true;
1985 if (current == Well::InjectorCMode::BHP) {
1986 bhp_controlled_well =
true;
1989 const Well::ProducerCMode& current = ws.production_cmode;
1990 if (current == Well::ProducerCMode::THP) {
1991 thp_controlled_well =
true;
1993 if (current == Well::ProducerCMode::BHP) {
1994 bhp_controlled_well =
true;
1997 if (!this->changed_to_open_this_step_ && (thp_controlled_well || bhp_controlled_well)) {
1999 double total_rate = 0.0;
2000 const double sign = this->isInjector() ? 1.0:-1.0;
2001 for (
int phase = 0; phase < np; ++phase){
2002 total_rate += sign * ws.surface_rates[phase];
2007 if (total_rate > 0) {
2008 for (
int phase = 0; phase < np; ++phase){
2009 well_potentials[phase] = sign * ws.surface_rates[phase];
2016 const auto& summaryState = ebosSimulator.vanguard().summaryState();
2017 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
2019 double bhp = this->mostStrictBhpFromBhpLimits(summaryState);
2027 if (this->isInjector())
2028 bhp = std::max(ws.bhp, bhp);
2030 bhp = std::min(ws.bhp, bhp);
2032 assert(std::abs(bhp) != std::numeric_limits<double>::max());
2033 computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
2036 well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
2039 const double sign = this->isInjector() ? 1.0:-1.0;
2040 double total_potential = 0.0;
2041 for (
int phase = 0; phase < np; ++phase){
2042 well_potentials[phase] *= sign;
2043 total_potential += well_potentials[phase];
2045 if (total_potential < 0.0 && this->param_.check_well_operability_) {
2047 this->operability_status_.has_negative_potentials =
true;
2048 const std::string msg = std::string(
"well ") + this->name() + std::string(
": has negative potentials and is not operable");
2049 deferred_logger.warning(
"NEGATIVE_POTENTIALS_INOPERABLE", msg);
2057 template<
typename TypeTag>
2062 this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
2063 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
2066 if constexpr (Base::has_polymermw) {
2067 if (this->isInjector()) {
2068 const auto& ws = well_state.well(this->index_of_well_);
2069 const auto& perf_data = ws.perf_data;
2070 const auto& water_velocity = perf_data.water_velocity;
2071 const auto& skin_pressure = perf_data.skin_pressure;
2072 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2073 this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
2074 this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
2078 for (
double v : this->primary_variables_) {
2080 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after update from wellState well: " << this->name(), deferred_logger);
2087 template<
typename TypeTag>
2089 StandardWell<TypeTag>::
2090 getRefDensity()
const
2092 return this->perf_densities_[0];
2098 template<
typename TypeTag>
2100 StandardWell<TypeTag>::
2101 updateWaterMobilityWithPolymer(
const Simulator& ebos_simulator,
2103 std::vector<EvalWell>& mob,
2104 DeferredLogger& deferred_logger)
const
2106 const int cell_idx = this->well_cells_[perf];
2107 const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
2108 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
2112 if (this->isInjector()) {
2114 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
2115 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2116 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
2119 if (PolymerModule::hasPlyshlog()) {
2122 if (this->isInjector() && this->wpolymer() == 0.) {
2127 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
2128 const EvalWell& bhp = this->getBhp();
2130 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
2131 double perf_dis_gas_rate = 0.;
2132 double perf_vap_oil_rate = 0.;
2133 double perf_vap_wat_rate = 0.;
2134 double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
2135 const double Tw = this->well_index_[perf] * trans_mult;
2136 computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
2137 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
2139 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
2140 const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
2141 const auto& scaled_drainage_info =
2142 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
2143 const double swcr = scaled_drainage_info.Swcr;
2144 const EvalWell poro = this->extendEval(int_quant.porosity());
2145 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
2147 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
2148 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2149 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
2151 if (PolymerModule::hasShrate()) {
2154 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
2156 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
2157 int_quant.pvtRegionIndex(),
2160 mob[waterCompIdx] /= shear_factor;
2164 template<
typename TypeTag>
2166 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
2173 typename SparseMatrixAdapter::MatrixBlock tmpMat;
2174 Dune::DynamicMatrix<Scalar> tmp;
2175 for (
auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
2177 const auto row_index = colC.index();
2179 for (
auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
2181 detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp);
2182 detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
2183 jacobian.addToBlock( row_index, colB.index(), tmpMat );
2189 template <
typename TypeTag>
2191 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
2192 const BVector& weights,
2193 const int pressureVarIndex,
2194 const bool use_well_weights,
2195 const WellState& well_state)
const
2210 int bhp_var_index = Bhp;
2212 auto cell_weights = weights[0];
2214 assert(this->duneC_.M() == weights.size());
2215 const int welldof_ind = this->duneC_.M() + this->index_of_well_;
2217 if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
2219 for (
auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
2220 const auto row_ind = colC.index();
2221 const auto& bw = weights[row_ind];
2223 assert((*colC).M() == bw.size());
2224 for (
size_t i = 0; i < bw.size(); ++i) {
2225 matel += (*colC)[bhp_var_index][i] * bw[i];
2228 jacobian[row_ind][welldof_ind] = matel;
2233 cell_weights /= nperf;
2235 BVectorWell bweights(1);
2236 size_t blockSz = this->numWellEq_;
2237 bweights[0].resize(blockSz);
2239 double diagElem = 0;
2241 if ( use_well_weights ){
2247 rhs[0].resize(blockSz);
2248 rhs[0][bhp_var_index] = 1.0;
2249 DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0];
2250 DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
2251 for (
size_t i = 0; i < blockSz; ++i) {
2253 for (
size_t j = 0; j < blockSz; ++j) {
2254 bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
2256 abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
2258 assert( abs_max > 0.0 );
2259 for (
size_t i = 0; i < blockSz; ++i) {
2260 bweights[0][i] /= abs_max;
2262 diagElem = 1.0/abs_max;
2265 if( this->isPressureControlled(well_state) ){
2266 bweights[0][blockSz-1] = 1.0;
2269 for (
size_t i = 0; i < cell_weights.size(); ++i) {
2270 bweights[0][i] = cell_weights[i];
2272 bweights[0][blockSz-1] = 0.0;
2274 const auto& locmat = this->duneD_[0][0];
2275 for (
size_t i = 0; i < cell_weights.size(); ++i) {
2276 diagElem += locmat[i][bhp_var_index]*cell_weights[i];
2283 jacobian[welldof_ind][welldof_ind] = diagElem;
2285 if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
2286 for (
auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
2287 const auto col_index = colB.index();
2288 const auto& bw = bweights[0];
2290 for (
size_t i = 0; i < bw.size(); ++i) {
2291 matel += (*colB)[i][pressureVarIndex] * bw[i];
2293 jacobian[welldof_ind][col_index] = matel;
2300 template<
typename TypeTag>
2301 typename StandardWell<TypeTag>::EvalWell
2302 StandardWell<TypeTag>::
2303 pskinwater(
const double throughput,
2304 const EvalWell& water_velocity,
2305 DeferredLogger& deferred_logger)
const
2307 if constexpr (Base::has_polymermw) {
2308 const int water_table_id = this->well_ecl_.getPolymerProperties().m_skprwattable;
2309 if (water_table_id <= 0) {
2310 OPM_DEFLOG_THROW(std::runtime_error,
"Unused SKPRWAT table id used for well " << name(), deferred_logger);
2312 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
2313 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2315 EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
2316 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
2319 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2320 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2328 template<
typename TypeTag>
2329 typename StandardWell<TypeTag>::EvalWell
2330 StandardWell<TypeTag>::
2331 pskin(
const double throughput,
2332 const EvalWell& water_velocity,
2333 const EvalWell& poly_inj_conc,
2334 DeferredLogger& deferred_logger)
const
2336 if constexpr (Base::has_polymermw) {
2337 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
2338 const EvalWell water_velocity_abs = abs(water_velocity);
2339 if (poly_inj_conc == 0.) {
2340 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
2342 const int polymer_table_id = this->well_ecl_.getPolymerProperties().m_skprpolytable;
2343 if (polymer_table_id <= 0) {
2344 OPM_DEFLOG_THROW(std::runtime_error,
"Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
2346 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
2347 const double reference_concentration = skprpolytable.refConcentration;
2348 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2350 EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
2351 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2352 if (poly_inj_conc == reference_concentration) {
2353 return sign * pskin_poly;
2356 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2357 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2358 return sign * pskin;
2360 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2361 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2369 template<
typename TypeTag>
2370 typename StandardWell<TypeTag>::EvalWell
2371 StandardWell<TypeTag>::
2372 wpolymermw(
const double throughput,
2373 const EvalWell& water_velocity,
2374 DeferredLogger& deferred_logger)
const
2376 if constexpr (Base::has_polymermw) {
2377 const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
2378 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2379 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2380 EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
2381 if (this->wpolymer() == 0.) {
2382 return molecular_weight;
2384 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2385 return molecular_weight;
2387 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2388 "while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
2396 template<
typename TypeTag>
2398 StandardWell<TypeTag>::
2399 updateWaterThroughput(
const double dt, WellState &well_state)
const
2401 if constexpr (Base::has_polymermw) {
2402 if (this->isInjector()) {
2403 auto& ws = well_state.well(this->index_of_well_);
2404 auto& perf_water_throughput = ws.perf_data.water_throughput;
2405 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2406 const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
2408 if (perf_water_vel > 0.) {
2409 perf_water_throughput[perf] += perf_water_vel * dt;
2420 template<
typename TypeTag>
2422 StandardWell<TypeTag>::
2423 handleInjectivityRate(
const Simulator& ebosSimulator,
2425 std::vector<EvalWell>& cq_s)
const
2427 const int cell_idx = this->well_cells_[perf];
2428 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2429 const auto& fs = int_quants.fluidState();
2430 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2431 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2432 const int wat_vel_index = Bhp + 1 + perf;
2433 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2437 cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
2443 template<
typename TypeTag>
2445 StandardWell<TypeTag>::
2446 handleInjectivityEquations(
const Simulator& ebosSimulator,
2447 const WellState& well_state,
2449 const EvalWell& water_flux_s,
2450 DeferredLogger& deferred_logger)
2452 const int cell_idx = this->well_cells_[perf];
2453 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2454 const auto& fs = int_quants.fluidState();
2455 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2456 const EvalWell water_flux_r = water_flux_s / b_w;
2457 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2458 const EvalWell water_velocity = water_flux_r / area;
2459 const int wat_vel_index = Bhp + 1 + perf;
2462 const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
2463 this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
2465 const auto& ws = well_state.well(this->index_of_well_);
2466 const auto& perf_data = ws.perf_data;
2467 const auto& perf_water_throughput = perf_data.water_throughput;
2468 const double throughput = perf_water_throughput[perf];
2469 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2471 EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
2472 poly_conc.setValue(this->wpolymer());
2475 const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
2476 - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
2478 this->resWell_[0][pskin_index] = eq_pskin.value();
2479 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
2480 this->duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
2481 this->duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
2485 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
2486 this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
2494 template<
typename TypeTag>
2496 StandardWell<TypeTag>::
2497 checkConvergenceExtraEqs(
const std::vector<double>& res,
2498 ConvergenceReport& report)
const
2503 if constexpr (Base::has_polymermw) {
2504 this->checkConvergencePolyMW(res, report, this->param_.max_residual_allowed_);
2512 template<
typename TypeTag>
2514 StandardWell<TypeTag>::
2515 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2516 const IntensiveQuantities& int_quants,
2517 const WellState& well_state,
2519 std::vector<RateVector>& connectionRates,
2520 DeferredLogger& deferred_logger)
const
2523 EvalWell cq_s_polymw = cq_s_poly;
2524 if (this->isInjector()) {
2525 const int wat_vel_index = Bhp + 1 + perf;
2526 const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
2527 if (water_velocity > 0.) {
2528 const auto& ws = well_state.well(this->index_of_well_);
2529 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2530 const double throughput = perf_water_throughput[perf];
2531 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2532 cq_s_polymw *= molecular_weight;
2538 }
else if (this->isProducer()) {
2539 if (cq_s_polymw < 0.) {
2540 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2547 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2555 template<
typename TypeTag>
2556 std::optional<double>
2557 StandardWell<TypeTag>::
2558 computeBhpAtThpLimitProd(
const WellState& well_state,
2559 const Simulator& ebos_simulator,
2560 const SummaryState& summary_state,
2561 DeferredLogger& deferred_logger)
const
2563 return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
2566 this->getALQ(well_state));
2569 template<
typename TypeTag>
2570 std::optional<double>
2571 StandardWell<TypeTag>::
2572 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
2573 const SummaryState& summary_state,
2574 DeferredLogger& deferred_logger,
2575 double alq_value)
const
2578 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2584 std::vector<double> rates(3);
2585 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2586 this->adaptRatesForVFP(rates);
2590 double max_pressure = 0.0;
2591 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2592 const int cell_idx = this->well_cells_[perf];
2593 const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
2594 const auto& fs = int_quants.fluidState();
2595 double pressure_cell = this->getPerfCellPressure(fs).value();
2596 max_pressure = std::max(max_pressure, pressure_cell);
2598 auto bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
2603 auto v = frates(*bhpAtLimit);
2604 if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }))
2607 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2611 std::vector<double> rates(3);
2612 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
2613 this->adaptRatesForVFP(rates);
2617 bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(fratesIter,
2622 v = frates(*bhpAtLimit);
2623 if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }))
2627 return std::nullopt;
2632 template<
typename TypeTag>
2633 std::optional<double>
2634 StandardWell<TypeTag>::
2635 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
2636 const SummaryState& summary_state,
2637 DeferredLogger& deferred_logger)
const
2640 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2646 std::vector<double> rates(3);
2647 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2651 return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitInj(frates,
2660 template<
typename TypeTag>
2662 StandardWell<TypeTag>::
2663 iterateWellEqWithControl(
const Simulator& ebosSimulator,
2665 const Well::InjectionControls& inj_controls,
2666 const Well::ProductionControls& prod_controls,
2667 WellState& well_state,
2668 const GroupState& group_state,
2669 DeferredLogger& deferred_logger)
2671 const int max_iter = this->param_.max_inner_iter_wells_;
2674 bool relax_convergence =
false;
2675 this->regularize_ =
false;
2677 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2679 if (it > this->param_.strict_inner_iter_wells_) {
2680 relax_convergence =
true;
2681 this->regularize_ =
true;
2684 auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
2686 converged = report.converged();
2692 solveEqAndUpdateWellState(well_state, deferred_logger);
2699 initPrimaryVariablesEvaluation();
2700 }
while (it < max_iter);
2706 template<
typename TypeTag>
2713 std::vector<double> well_q_s(this->num_components_, 0.);
2714 const EvalWell& bhp = this->getBhp();
2715 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2716 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2717 const int cell_idx = this->well_cells_[perf];
2718 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2719 std::vector<Scalar> mob(this->num_components_, 0.);
2720 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
2721 std::vector<Scalar> cq_s(this->num_components_, 0.);
2722 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
2723 const double Tw = this->well_index_[perf] * trans_mult;
2724 computePerfRateScalar(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2725 cq_s, deferred_logger);
2726 for (
int comp = 0; comp < this->num_components_; ++comp) {
2727 well_q_s[comp] += cq_s[comp];
2730 const auto& comm = this->parallel_well_info_.communication();
2731 if (comm.size() > 1)
2733 comm.sum(well_q_s.data(), well_q_s.size());
2742 template <
typename TypeTag>
2746 const std::function<
double(
const double)>& connPICalc,
2747 const std::vector<EvalWell>& mobility,
2748 double* connPI)
const
2751 const int np = this->number_of_phases_;
2752 for (
int p = 0; p < np; ++p) {
2755 const auto connMob =
2756 mobility[ this->flowPhaseToEbosCompIdx(p) ].value()
2757 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
2759 connPI[p] = connPICalc(connMob);
2762 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2763 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2765 const auto io = pu.phase_pos[Oil];
2766 const auto ig = pu.phase_pos[Gas];
2768 const auto vapoil = connPI[ig] * fs.Rv().value();
2769 const auto disgas = connPI[io] * fs.Rs().value();
2771 connPI[io] += vapoil;
2772 connPI[ig] += disgas;
2780 template <
typename TypeTag>
2782 StandardWell<TypeTag>::
2783 computeConnLevelInjInd(
const typename StandardWell<TypeTag>::FluidState& fs,
2784 const Phase preferred_phase,
2785 const std::function<
double(
const double)>& connIICalc,
2786 const std::vector<EvalWell>& mobility,
2788 DeferredLogger& deferred_logger)
const
2794 if (preferred_phase == Phase::GAS) {
2795 phase_pos = pu.phase_pos[Gas];
2797 else if (preferred_phase == Phase::OIL) {
2798 phase_pos = pu.phase_pos[Oil];
2800 else if (preferred_phase == Phase::WATER) {
2801 phase_pos = pu.phase_pos[Water];
2804 OPM_DEFLOG_THROW(NotImplemented,
2805 "Unsupported Injector Type ("
2806 <<
static_cast<int>(preferred_phase)
2807 <<
") for well " << this->name()
2808 <<
" during connection I.I. calculation",
2812 const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
2813 const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
2814 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:63
Collect per-connection static information to enable calculating connection-level or well-level produc...
Definition: WellProdIndexCalculator.hpp:36
double connectionProdIndStandard(const std::size_t connIdx, const double connMobility) const
Compute connection-level steady-state productivity index value using dynamic phase mobility.
Definition: WellProdIndexCalculator.cpp:110
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:37