22 #include <opm/common/utility/numeric/RootFinders.hpp>
23 #include <opm/input/eclipse/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>
87 StandardWell<TypeTag>::
88 computePerfRateEval(
const IntensiveQuantities& intQuants,
89 const std::vector<EvalWell>& mob,
94 std::vector<EvalWell>& cq_s,
95 double& perf_dis_gas_rate,
96 double& perf_vap_oil_rate,
97 DeferredLogger& deferred_logger)
const
99 const auto& fs = intQuants.fluidState();
100 const EvalWell pressure = this->extendEval(this->getPerfCellPressure(fs));
101 const EvalWell rs = this->extendEval(fs.Rs());
102 const EvalWell rv = this->extendEval(fs.Rv());
103 std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
104 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
105 if (!FluidSystem::phaseIsActive(phaseIdx)) {
108 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
109 b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
111 if constexpr (has_solvent) {
112 b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
115 if constexpr (has_zFraction) {
116 if (this->isInjector()) {
117 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
118 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
119 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
123 EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
125 if (this->isInjector()) {
126 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
127 skin_pressure = this->primary_variables_evaluation_[pskin_index];
132 std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
133 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
134 cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
154 template<
typename TypeTag>
156 StandardWell<TypeTag>::
157 computePerfRateScalar(
const IntensiveQuantities& intQuants,
158 const std::vector<Scalar>& mob,
163 std::vector<Scalar>& cq_s,
164 DeferredLogger& deferred_logger)
const
166 const auto& fs = intQuants.fluidState();
167 const Scalar pressure = this->getPerfCellPressure(fs).value();
168 const Scalar rs = fs.Rs().value();
169 const Scalar rv = fs.Rv().value();
170 std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
171 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
172 if (!FluidSystem::phaseIsActive(phaseIdx)) {
175 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
176 b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
178 if constexpr (has_solvent) {
179 b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
182 if constexpr (has_zFraction) {
183 if (this->isInjector()) {
184 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
185 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
186 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
190 Scalar skin_pressure =0.0;
192 if (this->isInjector()) {
193 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
194 skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
198 Scalar perf_dis_gas_rate = 0.0;
199 Scalar perf_vap_oil_rate = 0.0;
202 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
203 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
204 cmix_s[componentIdx] = getValue(this->wellSurfaceVolumeFraction(componentIdx));
224 template<
typename TypeTag>
225 template<
class Value>
227 StandardWell<TypeTag>::
228 computePerfRate(
const std::vector<Value>& mob,
229 const Value& pressure,
233 std::vector<Value>& b_perfcells_dense,
237 const Value& skin_pressure,
238 const std::vector<Value>& cmix_s,
239 std::vector<Value>& cq_s,
240 double& perf_dis_gas_rate,
241 double& perf_vap_oil_rate,
242 DeferredLogger& deferred_logger)
const
245 const Value well_pressure = bhp + this->perf_pressure_diffs_[perf];
246 Value drawdown = pressure - well_pressure;
247 if (this->isInjector()) {
248 drawdown += skin_pressure;
252 if ( drawdown > 0 ) {
254 if (!allow_cf && this->isInjector()) {
259 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
260 const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
261 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
264 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
265 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
266 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
267 const Value cq_sOil = cq_s[oilCompIdx];
268 const Value cq_sGas = cq_s[gasCompIdx];
269 const Value dis_gas = rs * cq_sOil;
270 const Value vap_oil = rv * cq_sGas;
272 cq_s[gasCompIdx] += dis_gas;
273 cq_s[oilCompIdx] += vap_oil;
276 if (this->isProducer()) {
277 perf_dis_gas_rate = getValue(dis_gas);
278 perf_vap_oil_rate = getValue(vap_oil);
284 if (!allow_cf && this->isProducer()) {
289 Value total_mob_dense = mob[0];
290 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
291 total_mob_dense += mob[componentIdx];
295 const Value cqt_i = - Tw * (total_mob_dense * drawdown);
298 Value volumeRatio = bhp * 0.0;
300 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
301 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
302 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
305 if constexpr (Indices::enableSolvent) {
306 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
309 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
310 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
311 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
313 const Value d = 1.0 - rv * rs;
316 std::ostringstream sstr;
317 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
318 <<
" during computePerfRate calculations with rs " << rs
319 <<
", rv " << rv <<
" and pressure " << pressure
320 <<
" obtaining d " << d
321 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
322 <<
" for this connection.";
323 deferred_logger.debug(sstr.str());
325 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
326 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
328 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
329 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
332 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
333 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
334 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
336 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
337 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
338 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
343 Value cqt_is = cqt_i/volumeRatio;
344 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
345 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
349 if (this->isProducer()) {
350 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
351 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
352 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
361 const double d = 1.0 - getValue(rv) * getValue(rs);
364 std::ostringstream sstr;
365 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
366 <<
" during computePerfRate calculations with rs " << rs
367 <<
", rv " << rv <<
" and pressure " << pressure
368 <<
" obtaining d " << d
369 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
370 <<
" for this connection.";
371 deferred_logger.debug(sstr.str());
375 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
378 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
386 template<
typename TypeTag>
388 StandardWell<TypeTag>::
389 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
391 const Well::InjectionControls& ,
392 const Well::ProductionControls& ,
393 WellState& well_state,
394 const GroupState& group_state,
395 DeferredLogger& deferred_logger)
399 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
404 this->invDuneD_ = 0.0;
405 this->resWell_ = 0.0;
407 assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
413 template<
typename TypeTag>
415 StandardWell<TypeTag>::
416 assembleWellEqWithoutIterationImpl(
const Simulator& ebosSimulator,
418 WellState& well_state,
419 const GroupState& group_state,
420 DeferredLogger& deferred_logger)
424 const double volume = 0.002831684659200;
425 auto& ws = well_state.well(this->index_of_well_);
427 ws.vaporized_oil_rate = 0;
428 ws.dissolved_gas_rate = 0;
430 const int np = this->number_of_phases_;
432 std::vector<RateVector> connectionRates = this->connectionRates_;
433 auto& perf_data = ws.perf_data;
434 auto& perf_rates = perf_data.phase_rates;
435 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
437 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
438 EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
439 EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
440 calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
443 if constexpr (has_polymer && Base::has_polymermw) {
444 if (this->isInjector()) {
445 handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
448 const int cell_idx = this->well_cells_[perf];
449 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
451 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
453 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
456 this->resWell_[0][componentIdx] += cq_s_effective.value();
459 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
461 this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq);
462 this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
465 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
466 this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
470 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
471 auto& perf_rate_solvent = perf_data.solvent_rates;
472 perf_rate_solvent[perf] = cq_s[componentIdx].value();
474 perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
478 if constexpr (has_zFraction) {
479 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
480 this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
485 this->connectionRates_ = connectionRates;
490 const auto& comm = this->parallel_well_info_.communication();
491 ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
492 ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
496 wellhelpers::sumDistributedWellEntries(this->invDuneD_[0][0], this->resWell_[0],
497 this->parallel_well_info_.communication());
499 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
502 EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
503 if (FluidSystem::numActivePhases() > 1) {
505 resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
507 resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
508 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
509 this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
511 this->resWell_[0][componentIdx] += resWell_loc.value();
514 const auto& summaryState = ebosSimulator.vanguard().summaryState();
515 const Schedule& schedule = ebosSimulator.vanguard().schedule();
516 this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger);
521 Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]);
523 OPM_DEFLOG_THROW(NumericalIssue,
"Error when inverting local well equations for well " + name(), deferred_logger);
530 template<
typename TypeTag>
532 StandardWell<TypeTag>::
533 calculateSinglePerf(
const Simulator& ebosSimulator,
535 WellState& well_state,
536 std::vector<RateVector>& connectionRates,
537 std::vector<EvalWell>& cq_s,
538 EvalWell& water_flux_s,
539 EvalWell& cq_s_zfrac_effective,
540 DeferredLogger& deferred_logger)
const
542 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
543 const EvalWell& bhp = this->getBhp();
544 const int cell_idx = this->well_cells_[perf];
545 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
546 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
547 getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
549 double perf_dis_gas_rate = 0.;
550 double perf_vap_oil_rate = 0.;
551 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
552 const double Tw = this->well_index_[perf] * trans_mult;
553 computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
554 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
556 auto& ws = well_state.well(this->index_of_well_);
557 auto& perf_data = ws.perf_data;
558 if constexpr (has_polymer && Base::has_polymermw) {
559 if (this->isInjector()) {
562 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
563 water_flux_s = cq_s[water_comp_idx];
566 handleInjectivityRate(ebosSimulator, perf, cq_s);
571 if (this->isProducer()) {
572 ws.dissolved_gas_rate += perf_dis_gas_rate;
573 ws.vaporized_oil_rate += perf_vap_oil_rate;
576 if constexpr (has_energy) {
577 connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
580 if constexpr (has_energy) {
582 auto fs = intQuants.fluidState();
583 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
584 if (!FluidSystem::phaseIsActive(phaseIdx)) {
589 EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
590 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
591 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
592 if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
593 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
596 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
597 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
602 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
604 std::ostringstream sstr;
605 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
606 <<
" during calculateSinglePerf with rs " << fs.Rs()
607 <<
", rv " << fs.Rv()
608 <<
" obtaining d " << d
609 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
610 <<
" for this connection.";
611 deferred_logger.debug(sstr.str());
612 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
614 if(FluidSystem::gasPhaseIdx == phaseIdx) {
615 cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
616 }
else if(FluidSystem::oilPhaseIdx == phaseIdx) {
618 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
624 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
626 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
627 fs.setTemperature(this->well_ecl_.temperature());
628 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
629 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
630 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
631 paramCache.setRegionIndex(pvtRegionIdx);
632 paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
633 paramCache.updatePhase(fs, phaseIdx);
635 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
636 fs.setDensity(phaseIdx, rho);
637 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
638 fs.setEnthalpy(phaseIdx, h);
641 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
642 connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
646 if constexpr (has_polymer) {
648 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
649 EvalWell cq_s_poly = cq_s[waterCompIdx];
650 if (this->isInjector()) {
651 cq_s_poly *= this->wpolymer();
653 cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
656 auto& perf_rate_polymer = perf_data.polymer_rates;
657 perf_rate_polymer[perf] = cq_s_poly.value();
659 cq_s_poly *= this->well_efficiency_factor_;
660 connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
662 if constexpr (Base::has_polymermw) {
663 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
667 if constexpr (has_foam) {
669 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
670 EvalWell cq_s_foam = cq_s[gasCompIdx] * this->well_efficiency_factor_;
671 if (this->isInjector()) {
672 cq_s_foam *= this->wfoam();
674 cq_s_foam *= this->extendEval(intQuants.foamConcentration());
676 connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
679 if constexpr (has_zFraction) {
681 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
682 cq_s_zfrac_effective = cq_s[gasCompIdx];
683 if (this->isInjector()) {
684 cq_s_zfrac_effective *= this->wsolvent();
685 }
else if (cq_s_zfrac_effective.value() != 0.0) {
686 const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
687 cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
689 auto& perf_rate_solvent = perf_data.solvent_rates;
690 perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
692 cq_s_zfrac_effective *= this->well_efficiency_factor_;
693 connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
696 if constexpr (has_brine) {
698 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
699 EvalWell cq_s_sm = cq_s[waterCompIdx];
700 if (this->isInjector()) {
701 cq_s_sm *= this->wsalt();
703 cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration());
706 auto& perf_rate_brine = perf_data.brine_rates;
707 perf_rate_brine[perf] = cq_s_sm.value();
709 cq_s_sm *= this->well_efficiency_factor_;
710 connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
713 if constexpr (has_micp) {
714 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
715 EvalWell cq_s_microbe = cq_s[waterCompIdx];
716 if (this->isInjector()) {
717 cq_s_microbe *= this->wmicrobes();
719 cq_s_microbe *= this->extendEval(intQuants.microbialConcentration());
721 connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe);
722 EvalWell cq_s_oxygen = cq_s[waterCompIdx];
723 if (this->isInjector()) {
724 cq_s_oxygen *= this->woxygen();
726 cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration());
728 connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen);
729 EvalWell cq_s_urea = cq_s[waterCompIdx];
730 if (this->isInjector()) {
731 cq_s_urea *= this->wurea();
733 cq_s_urea *= this->extendEval(intQuants.ureaConcentration());
735 connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea);
739 perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf];
745 template<
typename TypeTag>
747 StandardWell<TypeTag>::
748 getMobilityEval(
const Simulator& ebosSimulator,
750 std::vector<EvalWell>& mob,
751 DeferredLogger& deferred_logger)
const
753 const int cell_idx = this->well_cells_[perf];
754 assert (
int(mob.size()) == this->num_components_);
755 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
756 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
760 const int satid = this->saturation_table_number_[perf] - 1;
761 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
762 if( satid == satid_elem ) {
764 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
765 if (!FluidSystem::phaseIsActive(phaseIdx)) {
769 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
770 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
773 mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
777 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
778 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
779 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
782 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
785 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
786 if (!FluidSystem::phaseIsActive(phaseIdx)) {
790 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
791 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
795 if constexpr (has_solvent) {
796 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
801 if constexpr (has_polymer) {
802 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
803 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
808 if constexpr (!Base::has_polymermw) {
809 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
814 template<
typename TypeTag>
816 StandardWell<TypeTag>::
817 getMobilityScalar(
const Simulator& ebosSimulator,
819 std::vector<Scalar>& mob,
820 DeferredLogger& deferred_logger)
const
822 const int cell_idx = this->well_cells_[perf];
823 assert (
int(mob.size()) == this->num_components_);
824 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
825 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
829 const int satid = this->saturation_table_number_[perf] - 1;
830 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
831 if( satid == satid_elem ) {
833 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
834 if (!FluidSystem::phaseIsActive(phaseIdx)) {
838 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
839 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
842 mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
846 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
847 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
848 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
851 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
854 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
855 if (!FluidSystem::phaseIsActive(phaseIdx)) {
859 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
860 mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx));
864 if constexpr (has_solvent) {
865 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
870 if constexpr (has_polymer) {
871 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
872 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
877 if constexpr (!Base::has_polymermw) {
878 std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
879 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
880 for (
size_t i = 0; i < mob.size(); ++i) {
881 mob[i] = getValue(mob_eval[i]);
889 template<
typename TypeTag>
891 StandardWell<TypeTag>::
892 updateWellState(
const BVectorWell& dwells,
893 WellState& well_state,
894 DeferredLogger& deferred_logger)
const
896 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
898 updatePrimaryVariablesNewton(dwells, well_state, deferred_logger);
900 updateWellStateFromPrimaryVariables(well_state, deferred_logger);
901 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
908 template<
typename TypeTag>
910 StandardWell<TypeTag>::
911 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
913 DeferredLogger& deferred_logger)
const
915 const double dFLimit = this->param_.dwell_fraction_max_;
916 const double dBHPLimit = this->param_.dbhp_max_rel_;
917 this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
919 updateExtraPrimaryVariables(dwells);
921 for (
double v : this->primary_variables_) {
923 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after newton update well: " << this->name(), deferred_logger);
932 template<
typename TypeTag>
934 StandardWell<TypeTag>::
935 updateExtraPrimaryVariables(
const BVectorWell& dwells)
const
938 if constexpr (Base::has_polymermw) {
939 this->updatePrimaryVariablesPolyMW(dwells);
947 template<
typename TypeTag>
949 StandardWell<TypeTag>::
950 updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger)
const
952 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger);
955 if constexpr (Base::has_polymermw) {
956 this->updateWellStateFromPrimaryVariablesPolyMW(well_state);
964 template<
typename TypeTag>
966 StandardWell<TypeTag>::
967 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
972 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
973 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
975 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
976 std::vector<Scalar> mob(this->num_components_, 0.0);
977 getMobilityScalar(ebos_simulator, perf, mob, deferred_logger);
979 const int cell_idx = this->well_cells_[perf];
980 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
981 const auto& fs = int_quantities.fluidState();
983 double p_r = this->getPerfCellPressure(fs).value();
986 std::vector<double> b_perf(this->num_components_);
987 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
988 if (!FluidSystem::phaseIsActive(phase)) {
991 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
992 b_perf[comp_idx] = fs.invB(phase).value();
996 const double h_perf = this->perf_pressure_diffs_[perf];
997 const double pressure_diff = p_r - h_perf;
1002 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1003 deferred_logger.debug(
"CROSSFLOW_IPR",
1004 "cross flow found when updateIPR for well " + name()
1005 +
" . The connection is ignored in IPR calculations");
1011 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1016 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1017 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1018 for (
int p = 0; p < this->number_of_phases_; ++p) {
1019 const double tw_mob = tw_perf * mob[p] * b_perf[p];
1020 ipr_a_perf[p] += tw_mob * pressure_diff;
1021 ipr_b_perf[p] += tw_mob;
1025 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1026 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1027 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1028 const double rs = (fs.Rs()).value();
1029 const double rv = (fs.Rv()).value();
1031 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1032 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1034 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1035 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1037 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1038 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1040 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1041 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1044 for (
int p = 0; p < this->number_of_phases_; ++p) {
1046 this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
1047 this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
1050 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1051 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1055 template<
typename TypeTag>
1057 StandardWell<TypeTag>::
1058 checkOperabilityUnderBHPLimit(
const WellState& well_state,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1060 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1061 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1064 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1065 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1069 for (
int p = 0; p < this->number_of_phases_; ++p) {
1070 const double ipr_rate = this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
1071 if ( (this->isProducer() && ipr_rate < 0.) || (this->isInjector() && ipr_rate > 0.) ) {
1072 this->operability_status_.operable_under_only_bhp_limit =
false;
1078 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1082 std::vector<double> well_rates_bhp_limit;
1083 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1085 this->adaptRatesForVFP(well_rates_bhp_limit);
1086 const double thp = this->calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
1087 const double thp_limit = this->getTHPConstraint(summaryState);
1088 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1089 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1100 this->operability_status_.operable_under_only_bhp_limit =
true;
1101 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1109 template<
typename TypeTag>
1111 StandardWell<TypeTag>::
1112 checkOperabilityUnderTHPLimit(
const Simulator& ebos_simulator,
const WellState& well_state, DeferredLogger& deferred_logger)
1114 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1115 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger)
1116 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1119 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1121 const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1122 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1124 const double thp_limit = this->getTHPConstraint(summaryState);
1125 if (this->isProducer() && *obtain_bhp < thp_limit) {
1126 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1127 +
" bars is SMALLER than thp limit "
1128 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1129 +
" bars as a producer for well " + name();
1130 deferred_logger.debug(msg);
1132 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1133 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1134 +
" bars is LARGER than thp limit "
1135 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1136 +
" bars as a injector for well " + name();
1137 deferred_logger.debug(msg);
1140 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1141 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1142 if (!this->wellIsStopped()) {
1143 const double thp_limit = this->getTHPConstraint(summaryState);
1144 deferred_logger.debug(
" could not find bhp value at thp limit "
1145 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1146 +
" bar for well " + name() +
", the well might need to be closed ");
1155 template<
typename TypeTag>
1157 StandardWell<TypeTag>::
1158 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1160 bool all_drawdown_wrong_direction =
true;
1162 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1163 const int cell_idx = this->well_cells_[perf];
1164 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1165 const auto& fs = intQuants.fluidState();
1167 const double pressure = this->getPerfCellPressure(fs).value();
1168 const double bhp = this->getBhp().value();
1171 const double well_pressure = bhp + this->perf_pressure_diffs_[perf];
1172 const double drawdown = pressure - well_pressure;
1177 if ( (drawdown < 0. && this->isInjector()) ||
1178 (drawdown > 0. && this->isProducer()) ) {
1179 all_drawdown_wrong_direction =
false;
1184 const auto& comm = this->parallel_well_info_.communication();
1185 if (comm.size() > 1)
1187 all_drawdown_wrong_direction =
1188 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1191 return all_drawdown_wrong_direction;
1197 template<
typename TypeTag>
1199 StandardWell<TypeTag>::
1200 canProduceInjectWithCurrentBhp(
const Simulator& ebos_simulator,
1201 const WellState& well_state,
1202 DeferredLogger& deferred_logger)
1204 const double bhp = well_state.well(this->index_of_well_).bhp;
1205 std::vector<double> well_rates;
1206 computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
1208 const double sign = (this->isProducer()) ? -1. : 1.;
1209 const double threshold = sign * std::numeric_limits<double>::min();
1211 bool can_produce_inject =
false;
1212 for (
const auto value : well_rates) {
1213 if (this->isProducer() && value < threshold) {
1214 can_produce_inject =
true;
1216 }
else if (this->isInjector() && value > threshold) {
1217 can_produce_inject =
true;
1222 if (!can_produce_inject) {
1223 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1226 return can_produce_inject;
1233 template<
typename TypeTag>
1235 StandardWell<TypeTag>::
1236 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1238 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1244 template<
typename TypeTag>
1246 StandardWell<TypeTag>::
1247 computePropertiesForWellConnectionPressures(
const Simulator& ebosSimulator,
1248 const WellState& well_state,
1249 std::vector<double>& b_perf,
1250 std::vector<double>& rsmax_perf,
1251 std::vector<double>& rvmax_perf,
1252 std::vector<double>& surf_dens_perf)
const
1254 const int nperf = this->number_of_perforations_;
1256 b_perf.resize(nperf * this->num_components_);
1257 surf_dens_perf.resize(nperf * this->num_components_);
1258 const auto& ws = well_state.well(this->index_of_well_);
1260 const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1261 const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1262 const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1265 if (oilPresent && gasPresent) {
1266 rsmax_perf.resize(nperf);
1267 rvmax_perf.resize(nperf);
1271 const auto& perf_press = ws.perf_data.pressure;
1272 auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp,
1276 for (
int perf = 0; perf < nperf; ++perf) {
1277 const int cell_idx = this->well_cells_[perf];
1278 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1279 const auto& fs = intQuants.fluidState();
1281 const double p_avg = (perf_press[perf] + p_above[perf])/2;
1282 const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1283 const double saltConcentration = fs.saltConcentration().value();
1286 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1287 b_perf[ waterCompIdx + perf * this->num_components_] =
1288 FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration);
1292 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1293 const int gaspos = gasCompIdx + perf * this->num_components_;
1296 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1297 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1299 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1302 rv = oilrate / gasrate;
1304 rv = std::min(rv, rvmax_perf[perf]);
1306 b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
1309 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1313 b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1318 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1319 const int oilpos = oilCompIdx + perf * this->num_components_;
1321 rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1322 const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1324 const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1327 rs = gasrate / oilrate;
1329 rs = std::min(rs, rsmax_perf[perf]);
1330 b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1332 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1335 b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1340 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1341 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1345 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1346 surf_dens_perf[this->num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() );
1350 if constexpr (has_solvent) {
1351 b_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1352 surf_dens_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity();
1361 template<
typename TypeTag>
1365 const std::vector<double>& B_avg,
1367 const bool relax_tolerance)
const
1371 assert((
int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1373 std::vector<double> res;
1376 this->param_.max_residual_allowed_,
1377 this->param_.tolerance_wells_,
1378 this->param_.relaxed_tolerance_flow_well_,
1382 checkConvergenceExtraEqs(res, report);
1391 template<
typename TypeTag>
1399 auto fluidState = [&ebosSimulator,
this](
const int perf)
1401 const auto cell_idx = this->well_cells_[perf];
1402 return ebosSimulator.model()
1403 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
1406 const int np = this->number_of_phases_;
1407 auto setToZero = [np](
double* x) ->
void
1409 std::fill_n(x, np, 0.0);
1412 auto addVector = [np](
const double* src,
double* dest) ->
void
1414 std::transform(src, src + np, dest, dest, std::plus<>{});
1417 auto& ws = well_state.well(this->index_of_well_);
1418 auto& perf_data = ws.perf_data;
1419 auto* wellPI = ws.productivity_index.data();
1420 auto* connPI = perf_data.prod_index.data();
1424 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1425 auto subsetPerfID = 0;
1427 for (
const auto& perf : *this->perf_data_) {
1428 auto allPerfID = perf.ecl_index;
1430 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1435 std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
1436 getMobilityEval(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1438 const auto& fs = fluidState(subsetPerfID);
1441 if (this->isInjector()) {
1442 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1443 mob, connPI, deferred_logger);
1446 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1449 addVector(connPI, wellPI);
1456 const auto& comm = this->parallel_well_info_.communication();
1457 if (comm.size() > 1) {
1458 comm.sum(wellPI, np);
1461 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1462 "Internal logic error in processing connections for PI/II");
1467 template<
typename TypeTag>
1469 StandardWell<TypeTag>::
1470 computeWellConnectionDensitesPressures(
const Simulator& ebosSimulator,
1471 const WellState& well_state,
1472 const std::vector<double>& b_perf,
1473 const std::vector<double>& rsmax_perf,
1474 const std::vector<double>& rvmax_perf,
1475 const std::vector<double>& surf_dens_perf,
1476 DeferredLogger& deferred_logger)
1479 const int nperf = this->number_of_perforations_;
1480 const int np = this->number_of_phases_;
1481 std::vector<double> perfRates(b_perf.size(),0.0);
1482 const auto& ws = well_state.well(this->index_of_well_);
1483 const auto& perf_data = ws.perf_data;
1484 const auto& perf_rates_state = perf_data.phase_rates;
1486 for (
int perf = 0; perf < nperf; ++perf) {
1487 for (
int comp = 0; comp < np; ++comp) {
1488 perfRates[perf * this->num_components_ + comp] = perf_rates_state[perf * np + this->ebosCompIdxToFlowCompIdx(comp)];
1492 if constexpr (has_solvent) {
1493 const auto& solvent_perf_rates_state = perf_data.solvent_rates;
1494 for (
int perf = 0; perf < nperf; ++perf) {
1495 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = solvent_perf_rates_state[perf];
1502 bool all_zero = std::all_of(perfRates.begin(), perfRates.end(), [](
double val) { return val == 0.0; });
1503 if ( all_zero && this->isProducer() ) {
1504 double total_tw = 0;
1505 for (
int perf = 0; perf < nperf; ++perf) {
1506 total_tw += this->well_index_[perf];
1508 for (
int perf = 0; perf < nperf; ++perf) {
1509 const int cell_idx = this->well_cells_[perf];
1510 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1511 const auto& fs = intQuants.fluidState();
1512 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1513 double total_mobility = 0.0;
1514 for (
int p = 0; p < np; ++p) {
1515 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1516 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1518 if constexpr (has_solvent) {
1519 total_mobility += intQuants.solventInverseFormationVolumeFactor().value() * intQuants.solventMobility().value();
1521 for (
int p = 0; p < np; ++p) {
1522 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1523 perfRates[perf * this->num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1525 if constexpr (has_solvent) {
1526 perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility;
1531 this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
1533 this->computeConnectionPressureDelta();
1540 template<
typename TypeTag>
1542 StandardWell<TypeTag>::
1543 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1544 const WellState& well_state,
1545 DeferredLogger& deferred_logger)
1550 std::vector<double> b_perf;
1551 std::vector<double> rsmax_perf;
1552 std::vector<double> rvmax_perf;
1553 std::vector<double> surf_dens_perf;
1554 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1555 computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
1562 template<
typename TypeTag>
1564 StandardWell<TypeTag>::
1565 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
1567 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1571 BVectorWell dx_well(1);
1572 dx_well[0].resize(this->numWellEq_);
1573 this->invDuneD_.mv(this->resWell_, dx_well);
1575 updateWellState(dx_well, well_state, deferred_logger);
1582 template<
typename TypeTag>
1584 StandardWell<TypeTag>::
1585 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1586 const WellState& well_state,
1587 DeferredLogger& deferred_logger)
1589 updatePrimaryVariables(well_state, deferred_logger);
1590 initPrimaryVariablesEvaluation();
1591 computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
1592 this->computeAccumWell();
1597 template<
typename TypeTag>
1600 apply(
const BVector& x, BVector& Ax)
const
1602 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1604 if (this->param_.matrix_add_well_contributions_)
1609 assert( this->Bx_.size() == this->duneB_.N() );
1610 assert( this->invDrw_.size() == this->invDuneD_.N() );
1613 this->parallelB_.mv(x, this->Bx_);
1618 BVectorWell& invDBx = this->invDrw_;
1619 this->invDuneD_.mv(this->Bx_, invDBx);
1622 this->duneC_.mmtv(invDBx,Ax);
1628 template<
typename TypeTag>
1631 apply(BVector& r)
const
1633 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1635 assert( this->invDrw_.size() == this->invDuneD_.N() );
1638 this->invDuneD_.mv(this->resWell_, this->invDrw_);
1640 this->duneC_.mmtv(this->invDrw_, r);
1643 template<
typename TypeTag>
1648 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1650 BVectorWell resWell = this->resWell_;
1652 this->parallelB_.mmv(x, resWell);
1654 this->invDuneD_.mv(resWell, xw);
1661 template<
typename TypeTag>
1668 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1671 xw[0].resize(this->numWellEq_);
1673 recoverSolutionWell(x, xw);
1674 updateWellState(xw, well_state, deferred_logger);
1680 template<
typename TypeTag>
1685 std::vector<double>& well_flux,
1689 const int np = this->number_of_phases_;
1690 well_flux.resize(np, 0.0);
1692 const bool allow_cf = this->getAllowCrossFlow();
1694 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1695 const int cell_idx = this->well_cells_[perf];
1696 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1698 std::vector<Scalar> mob(this->num_components_, 0.);
1699 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
1700 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
1701 const double Tw = this->well_index_[perf] * trans_mult;
1703 std::vector<Scalar> cq_s(this->num_components_, 0.);
1704 computePerfRateScalar(intQuants, mob, bhp, Tw, perf, allow_cf,
1705 cq_s, deferred_logger);
1707 for(
int p = 0; p < np; ++p) {
1708 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
1711 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1716 template<
typename TypeTag>
1718 StandardWell<TypeTag>::
1719 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
1721 std::vector<double>& well_flux,
1722 DeferredLogger& deferred_logger)
const
1728 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1729 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1730 auto& ws = well_state_copy.well(this->index_of_well_);
1733 if (this->well_ecl_.isInjector()) {
1734 ws.injection_cmode = Well::InjectorCMode::BHP;
1736 ws.production_cmode = Well::ProducerCMode::BHP;
1741 const int np = this->number_of_phases_;
1742 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1743 for (
int phase = 0; phase < np; ++phase){
1744 well_state_copy.wellRates(this->index_of_well_)[phase]
1745 = sign * ws.well_potentials[phase];
1749 StandardWell<TypeTag> well(*
this);
1750 well.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1752 const double dt = ebosSimulator.timeStepSize();
1753 bool converged = well.iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
1755 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1756 " potentials are computed based on unconverged solution";
1757 deferred_logger.debug(msg);
1759 well.updatePrimaryVariables(well_state_copy, deferred_logger);
1760 well.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
1761 well.initPrimaryVariablesEvaluation();
1762 well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1768 template<
typename TypeTag>
1770 StandardWell<TypeTag>::
1771 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
1772 DeferredLogger& deferred_logger,
1773 const WellState &well_state)
const
1775 std::vector<double> potentials(this->number_of_phases_, 0.0);
1776 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1778 const auto& well = this->well_ecl_;
1779 if (well.isInjector()){
1780 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1781 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1782 if (bhp_at_thp_limit) {
1783 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1784 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1786 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1787 "Failed in getting converged thp based potential calculation for well "
1788 + name() +
". Instead the bhp based value is used");
1789 const double bhp = controls.bhp_limit;
1790 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1793 computeWellRatesWithThpAlqProd(
1794 ebos_simulator, summary_state,
1795 deferred_logger, potentials, this->getALQ(well_state)
1802 template<
typename TypeTag>
1804 StandardWell<TypeTag>::
1805 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &ebos_simulator,
1806 const SummaryState &summary_state,
1807 DeferredLogger &deferred_logger,
1808 std::vector<double> &potentials,
1812 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1813 ebos_simulator, summary_state, deferred_logger, alq);
1814 if (bhp_at_thp_limit) {
1815 const auto& controls = this->well_ecl_.productionControls(summary_state);
1816 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1817 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1820 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1821 "Failed in getting converged thp based potential calculation for well "
1822 + name() +
". Instead the bhp based value is used");
1823 const auto& controls = this->well_ecl_.productionControls(summary_state);
1824 bhp = controls.bhp_limit;
1825 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1830 template<
typename TypeTag>
1832 StandardWell<TypeTag>::
1833 computeWellRatesWithThpAlqProd(
const Simulator &ebos_simulator,
1834 const SummaryState &summary_state,
1835 DeferredLogger &deferred_logger,
1836 std::vector<double> &potentials,
1840 computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1847 template<
typename TypeTag>
1852 std::vector<double>& well_potentials,
1855 const int np = this->number_of_phases_;
1856 well_potentials.resize(np, 0.0);
1858 if (this->wellIsStopped()) {
1862 this->operability_status_.has_negative_potentials =
false;
1864 bool thp_controlled_well =
false;
1865 bool bhp_controlled_well =
false;
1866 const auto& ws = well_state.well(this->index_of_well_);
1867 if (this->isInjector()) {
1868 const Well::InjectorCMode& current = ws.injection_cmode;
1869 if (current == Well::InjectorCMode::THP) {
1870 thp_controlled_well =
true;
1872 if (current == Well::InjectorCMode::BHP) {
1873 bhp_controlled_well =
true;
1876 const Well::ProducerCMode& current = ws.production_cmode;
1877 if (current == Well::ProducerCMode::THP) {
1878 thp_controlled_well =
true;
1880 if (current == Well::ProducerCMode::BHP) {
1881 bhp_controlled_well =
true;
1884 if (thp_controlled_well || bhp_controlled_well) {
1886 double total_rate = 0.0;
1887 const double sign = this->isInjector() ? 1.0:-1.0;
1888 for (
int phase = 0; phase < np; ++phase){
1889 total_rate += sign * ws.surface_rates[phase];
1894 if (total_rate > 0) {
1895 for (
int phase = 0; phase < np; ++phase){
1896 well_potentials[phase] = sign * ws.surface_rates[phase];
1903 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1904 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1906 double bhp = this->mostStrictBhpFromBhpLimits(summaryState);
1914 if (this->isInjector())
1915 bhp = std::max(ws.bhp, bhp);
1917 bhp = std::min(ws.bhp, bhp);
1919 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1920 computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
1923 well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
1926 const double sign = this->isInjector() ? 1.0:-1.0;
1927 double total_potential = 0.0;
1928 for (
int phase = 0; phase < np; ++phase){
1929 well_potentials[phase] *= sign;
1930 total_potential += well_potentials[phase];
1932 if (total_potential < 0.0 && this->param_.check_well_operability_) {
1934 this->operability_status_.has_negative_potentials =
true;
1935 const std::string msg = std::string(
"well ") + this->name() + std::string(
": has negative potentials and is not operable");
1936 deferred_logger.warning(
"NEGATIVE_POTENTIALS_INOPERABLE", msg);
1944 template<
typename TypeTag>
1949 this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
1950 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1953 if constexpr (Base::has_polymermw) {
1954 if (this->isInjector()) {
1955 const auto& ws = well_state.well(this->index_of_well_);
1956 const auto& perf_data = ws.perf_data;
1957 const auto& water_velocity = perf_data.water_velocity;
1958 const auto& skin_pressure = perf_data.skin_pressure;
1959 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1960 this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
1961 this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
1965 for (
double v : this->primary_variables_) {
1967 OPM_DEFLOG_THROW(NumericalIssue,
"Infinite primary variable after update from wellState well: " << this->name(), deferred_logger);
1974 template<
typename TypeTag>
1976 StandardWell<TypeTag>::
1977 getRefDensity()
const
1979 return this->perf_densities_[0];
1985 template<
typename TypeTag>
1987 StandardWell<TypeTag>::
1988 updateWaterMobilityWithPolymer(
const Simulator& ebos_simulator,
1990 std::vector<EvalWell>& mob,
1991 DeferredLogger& deferred_logger)
const
1993 const int cell_idx = this->well_cells_[perf];
1994 const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1995 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1999 if (this->isInjector()) {
2001 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
2002 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2003 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
2006 if (PolymerModule::hasPlyshlog()) {
2009 if (this->isInjector() && this->wpolymer() == 0.) {
2014 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
2015 const EvalWell& bhp = this->getBhp();
2017 std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
2018 double perf_dis_gas_rate = 0.;
2019 double perf_vap_oil_rate = 0.;
2020 double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
2021 const double Tw = this->well_index_[perf] * trans_mult;
2022 computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
2023 cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
2025 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
2026 const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
2027 const auto& scaled_drainage_info =
2028 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
2029 const double swcr = scaled_drainage_info.Swcr;
2030 const EvalWell poro = this->extendEval(int_quant.porosity());
2031 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
2033 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
2034 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2035 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
2037 if (PolymerModule::hasShrate()) {
2040 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
2042 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
2043 int_quant.pvtRegionIndex(),
2046 mob[waterCompIdx] /= shear_factor;
2050 template<
typename TypeTag>
2052 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
2059 typename SparseMatrixAdapter::MatrixBlock tmpMat;
2060 Dune::DynamicMatrix<Scalar> tmp;
2061 for (
auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
2063 const auto row_index = colC.index();
2065 for (
auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
2067 Detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp);
2068 Detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
2069 jacobian.addToBlock( row_index, colB.index(), tmpMat );
2078 template<
typename TypeTag>
2079 typename StandardWell<TypeTag>::EvalWell
2080 StandardWell<TypeTag>::
2081 pskinwater(
const double throughput,
2082 const EvalWell& water_velocity,
2083 DeferredLogger& deferred_logger)
const
2085 if constexpr (Base::has_polymermw) {
2086 const int water_table_id = this->well_ecl_.getPolymerProperties().m_skprwattable;
2087 if (water_table_id <= 0) {
2088 OPM_DEFLOG_THROW(std::runtime_error,
"Unused SKPRWAT table id used for well " << name(), deferred_logger);
2090 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
2091 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2093 EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
2094 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
2097 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2098 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2106 template<
typename TypeTag>
2107 typename StandardWell<TypeTag>::EvalWell
2108 StandardWell<TypeTag>::
2109 pskin(
const double throughput,
2110 const EvalWell& water_velocity,
2111 const EvalWell& poly_inj_conc,
2112 DeferredLogger& deferred_logger)
const
2114 if constexpr (Base::has_polymermw) {
2115 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
2116 const EvalWell water_velocity_abs = abs(water_velocity);
2117 if (poly_inj_conc == 0.) {
2118 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
2120 const int polymer_table_id = this->well_ecl_.getPolymerProperties().m_skprpolytable;
2121 if (polymer_table_id <= 0) {
2122 OPM_DEFLOG_THROW(std::runtime_error,
"Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
2124 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
2125 const double reference_concentration = skprpolytable.refConcentration;
2126 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2128 EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
2129 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2130 if (poly_inj_conc == reference_concentration) {
2131 return sign * pskin_poly;
2134 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2135 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2136 return sign * pskin;
2138 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2139 "while injecting skin pressure is requested for well " << name(), deferred_logger);
2147 template<
typename TypeTag>
2148 typename StandardWell<TypeTag>::EvalWell
2149 StandardWell<TypeTag>::
2150 wpolymermw(
const double throughput,
2151 const EvalWell& water_velocity,
2152 DeferredLogger& deferred_logger)
const
2154 if constexpr (Base::has_polymermw) {
2155 const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
2156 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2157 const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2158 EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
2159 if (this->wpolymer() == 0.) {
2160 return molecular_weight;
2162 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2163 return molecular_weight;
2165 OPM_DEFLOG_THROW(std::runtime_error,
"Polymermw is not activated, "
2166 "while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
2174 template<
typename TypeTag>
2176 StandardWell<TypeTag>::
2177 updateWaterThroughput(
const double dt, WellState &well_state)
const
2179 if constexpr (Base::has_polymermw) {
2180 if (this->isInjector()) {
2181 auto& ws = well_state.well(this->index_of_well_);
2182 auto& perf_water_throughput = ws.perf_data.water_throughput;
2183 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2184 const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
2186 if (perf_water_vel > 0.) {
2187 perf_water_throughput[perf] += perf_water_vel * dt;
2198 template<
typename TypeTag>
2200 StandardWell<TypeTag>::
2201 handleInjectivityRate(
const Simulator& ebosSimulator,
2203 std::vector<EvalWell>& cq_s)
const
2205 const int cell_idx = this->well_cells_[perf];
2206 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2207 const auto& fs = int_quants.fluidState();
2208 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2209 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2210 const int wat_vel_index = Bhp + 1 + perf;
2211 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2215 cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
2221 template<
typename TypeTag>
2223 StandardWell<TypeTag>::
2224 handleInjectivityEquations(
const Simulator& ebosSimulator,
2225 const WellState& well_state,
2227 const EvalWell& water_flux_s,
2228 DeferredLogger& deferred_logger)
2230 const int cell_idx = this->well_cells_[perf];
2231 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2232 const auto& fs = int_quants.fluidState();
2233 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2234 const EvalWell water_flux_r = water_flux_s / b_w;
2235 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2236 const EvalWell water_velocity = water_flux_r / area;
2237 const int wat_vel_index = Bhp + 1 + perf;
2240 const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
2241 this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
2243 const auto& ws = well_state.well(this->index_of_well_);
2244 const auto& perf_data = ws.perf_data;
2245 const auto& perf_water_throughput = perf_data.water_throughput;
2246 const double throughput = perf_water_throughput[perf];
2247 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2249 EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
2250 poly_conc.setValue(this->wpolymer());
2253 const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
2254 - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
2256 this->resWell_[0][pskin_index] = eq_pskin.value();
2257 for (
int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
2258 this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
2259 this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
2263 for (
int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
2264 this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
2272 template<
typename TypeTag>
2274 StandardWell<TypeTag>::
2275 checkConvergenceExtraEqs(
const std::vector<double>& res,
2276 ConvergenceReport& report)
const
2281 if constexpr (Base::has_polymermw) {
2282 this->checkConvergencePolyMW(res, report, this->param_.max_residual_allowed_);
2290 template<
typename TypeTag>
2292 StandardWell<TypeTag>::
2293 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2294 const IntensiveQuantities& int_quants,
2295 const WellState& well_state,
2297 std::vector<RateVector>& connectionRates,
2298 DeferredLogger& deferred_logger)
const
2301 EvalWell cq_s_polymw = cq_s_poly;
2302 if (this->isInjector()) {
2303 const int wat_vel_index = Bhp + 1 + perf;
2304 const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
2305 if (water_velocity > 0.) {
2306 const auto& ws = well_state.well(this->index_of_well_);
2307 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2308 const double throughput = perf_water_throughput[perf];
2309 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2310 cq_s_polymw *= molecular_weight;
2316 }
else if (this->isProducer()) {
2317 if (cq_s_polymw < 0.) {
2318 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2325 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2333 template<
typename TypeTag>
2334 std::optional<double>
2335 StandardWell<TypeTag>::
2336 computeBhpAtThpLimitProd(
const WellState& well_state,
2337 const Simulator& ebos_simulator,
2338 const SummaryState& summary_state,
2339 DeferredLogger& deferred_logger)
const
2341 return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
2344 this->getALQ(well_state));
2347 template<
typename TypeTag>
2348 std::optional<double>
2349 StandardWell<TypeTag>::
2350 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
2351 const SummaryState& summary_state,
2352 DeferredLogger& deferred_logger,
2353 double alq_value)
const
2356 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2362 std::vector<double> rates(3);
2363 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2364 this->adaptRatesForVFP(rates);
2368 double max_pressure = 0.0;
2369 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2370 const int cell_idx = this->well_cells_[perf];
2371 const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
2372 const auto& fs = int_quants.fluidState();
2373 double pressure_cell = this->getPerfCellPressure(fs).value();
2374 max_pressure = std::max(max_pressure, pressure_cell);
2376 auto bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
2381 auto v = frates(*bhpAtLimit);
2382 if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }))
2385 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2389 std::vector<double> rates(3);
2390 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
2394 bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(fratesIter,
2399 v = frates(*bhpAtLimit);
2400 if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }))
2404 return std::nullopt;
2409 template<
typename TypeTag>
2410 std::optional<double>
2411 StandardWell<TypeTag>::
2412 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
2413 const SummaryState& summary_state,
2414 DeferredLogger& deferred_logger)
const
2417 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2423 std::vector<double> rates(3);
2424 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2428 return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitInj(frates,
2437 template<
typename TypeTag>
2439 StandardWell<TypeTag>::
2440 iterateWellEqWithControl(
const Simulator& ebosSimulator,
2442 const Well::InjectionControls& inj_controls,
2443 const Well::ProductionControls& prod_controls,
2444 WellState& well_state,
2445 const GroupState& group_state,
2446 DeferredLogger& deferred_logger)
2448 const int max_iter = this->param_.max_inner_iter_wells_;
2452 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2454 auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger);
2456 converged = report.converged();
2462 solveEqAndUpdateWellState(well_state, deferred_logger);
2469 initPrimaryVariablesEvaluation();
2470 }
while (it < max_iter);
2476 template<
typename TypeTag>
2483 std::vector<double> well_q_s(this->num_components_, 0.);
2484 const EvalWell& bhp = this->getBhp();
2485 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2486 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2487 const int cell_idx = this->well_cells_[perf];
2488 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2489 std::vector<Scalar> mob(this->num_components_, 0.);
2490 getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
2491 std::vector<Scalar> cq_s(this->num_components_, 0.);
2492 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
2493 const double Tw = this->well_index_[perf] * trans_mult;
2494 computePerfRateScalar(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2495 cq_s, deferred_logger);
2496 for (
int comp = 0; comp < this->num_components_; ++comp) {
2497 well_q_s[comp] += cq_s[comp];
2500 const auto& comm = this->parallel_well_info_.communication();
2501 if (comm.size() > 1)
2503 comm.sum(well_q_s.data(), well_q_s.size());
2512 template <
typename TypeTag>
2516 const std::function<
double(
const double)>& connPICalc,
2517 const std::vector<EvalWell>& mobility,
2518 double* connPI)
const
2521 const int np = this->number_of_phases_;
2522 for (
int p = 0; p < np; ++p) {
2525 const auto connMob =
2526 mobility[ this->flowPhaseToEbosCompIdx(p) ].value()
2527 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
2529 connPI[p] = connPICalc(connMob);
2532 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2533 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2535 const auto io = pu.phase_pos[Oil];
2536 const auto ig = pu.phase_pos[Gas];
2538 const auto vapoil = connPI[ig] * fs.Rv().value();
2539 const auto disgas = connPI[io] * fs.Rs().value();
2541 connPI[io] += vapoil;
2542 connPI[ig] += disgas;
2550 template <
typename TypeTag>
2552 StandardWell<TypeTag>::
2553 computeConnLevelInjInd(
const typename StandardWell<TypeTag>::FluidState& fs,
2554 const Phase preferred_phase,
2555 const std::function<
double(
const double)>& connIICalc,
2556 const std::vector<EvalWell>& mobility,
2558 DeferredLogger& deferred_logger)
const
2564 if (preferred_phase == Phase::GAS) {
2565 phase_pos = pu.phase_pos[Gas];
2567 else if (preferred_phase == Phase::OIL) {
2568 phase_pos = pu.phase_pos[Oil];
2570 else if (preferred_phase == Phase::WATER) {
2571 phase_pos = pu.phase_pos[Water];
2574 OPM_DEFLOG_THROW(NotImplemented,
2575 "Unsupported Injector Type ("
2576 <<
static_cast<int>(preferred_phase)
2577 <<
") for well " << this->name()
2578 <<
" during connection I.I. calculation",
2582 const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
2583 const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
2584 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: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:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33