22#include <opm/common/Exceptions.hpp>
24#include <opm/input/eclipse/Units/Units.hpp>
26#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
27#include <opm/simulators/wells/StandardWellAssemble.hpp>
28#include <opm/simulators/wells/VFPHelpers.hpp>
29#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
30#include <opm/simulators/wells/WellConvergence.hpp>
32#include <fmt/format.h>
41 template<
typename TypeTag>
42 StandardWell<TypeTag>::
43 StandardWell(
const Well& well,
44 const ParallelWellInfo& pw_info,
46 const ModelParameters& param,
47 const RateConverterType& rate_converter,
48 const int pvtRegionIdx,
49 const int num_components,
51 const int index_of_well,
52 const std::vector<PerforationData>& perf_data)
53 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
54 , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
57 assert(this->num_components_ == numWellConservationEq);
64 template<
typename TypeTag>
66 StandardWell<TypeTag>::
67 init(
const PhaseUsage* phase_usage_arg,
68 const std::vector<double>& depth_arg,
69 const double gravity_arg,
71 const std::vector< Scalar >& B_avg,
72 const bool changed_to_open_this_step)
74 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
75 this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
82 template<
typename TypeTag>
83 void StandardWell<TypeTag>::
84 initPrimaryVariablesEvaluation()
86 this->primary_variables_.init();
93 template<
typename TypeTag>
95 StandardWell<TypeTag>::
96 computePerfRateEval(
const IntensiveQuantities& intQuants,
97 const std::vector<EvalWell>& mob,
102 std::vector<EvalWell>& cq_s,
103 double& perf_dis_gas_rate,
104 double& perf_dis_gas_rate_in_water,
105 double& perf_vap_oil_rate,
106 double& perf_vap_wat_rate,
107 DeferredLogger& deferred_logger)
const
109 const auto& fs = intQuants.fluidState();
110 const EvalWell pressure = this->extendEval(this->getPerfCellPressure(fs));
111 const EvalWell rs = this->extendEval(fs.Rs());
112 const EvalWell rv = this->extendEval(fs.Rv());
113 const EvalWell rvw = this->extendEval(fs.Rvw());
114 const EvalWell rsw = this->extendEval(fs.Rsw());
117 std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
118 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
119 if (!FluidSystem::phaseIsActive(phaseIdx)) {
122 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
123 b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
125 if constexpr (has_solvent) {
126 b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
129 if constexpr (has_zFraction) {
130 if (this->isInjector()) {
131 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
132 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
133 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
137 EvalWell skin_pressure = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
139 if (this->isInjector()) {
140 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
141 skin_pressure = this->primary_variables_.eval(pskin_index);
146 std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->primary_variables_.numWellEq() + Indices::numEq});
147 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
148 cmix_s[componentIdx] = this->primary_variables_.surfaceVolumeFraction(componentIdx);
166 perf_dis_gas_rate_in_water,
172 template<
typename TypeTag>
174 StandardWell<TypeTag>::
175 computePerfRateScalar(
const IntensiveQuantities& intQuants,
176 const std::vector<Scalar>& mob,
181 std::vector<Scalar>& cq_s,
182 DeferredLogger& deferred_logger)
const
184 const auto& fs = intQuants.fluidState();
185 const Scalar pressure = this->getPerfCellPressure(fs).value();
186 const Scalar rs = fs.Rs().value();
187 const Scalar rv = fs.Rv().value();
188 const Scalar rvw = fs.Rvw().value();
189 const Scalar rsw = fs.Rsw().value();
190 std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
191 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
192 if (!FluidSystem::phaseIsActive(phaseIdx)) {
195 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
196 b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
198 if constexpr (has_solvent) {
199 b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
202 if constexpr (has_zFraction) {
203 if (this->isInjector()) {
204 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
205 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
206 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
210 Scalar skin_pressure =0.0;
212 if (this->isInjector()) {
213 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
214 skin_pressure = getValue(this->primary_variables_.eval(pskin_index));
218 Scalar perf_dis_gas_rate = 0.0;
219 Scalar perf_vap_oil_rate = 0.0;
220 Scalar perf_vap_wat_rate = 0.0;
221 Scalar perf_dis_gas_rate_in_water = 0.0;
224 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
225 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
226 cmix_s[componentIdx] = getValue(this->primary_variables_.surfaceVolumeFraction(componentIdx));
244 perf_dis_gas_rate_in_water,
250 template<
typename TypeTag>
251 template<
class Value>
253 StandardWell<TypeTag>::
254 computePerfRate(
const std::vector<Value>& mob,
255 const Value& pressure,
261 std::vector<Value>& b_perfcells_dense,
265 const Value& skin_pressure,
266 const std::vector<Value>& cmix_s,
267 std::vector<Value>& cq_s,
268 double& perf_dis_gas_rate,
269 double& perf_dis_gas_rate_in_water,
270 double& perf_vap_oil_rate,
271 double& perf_vap_wat_rate,
272 DeferredLogger& deferred_logger)
const
275 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
276 Value drawdown = pressure - well_pressure;
277 if (this->isInjector()) {
278 drawdown += skin_pressure;
282 if ( drawdown > 0 ) {
284 if (!allow_cf && this->isInjector()) {
289 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
290 const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
291 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
294 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
295 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
296 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
297 const Value cq_sOil = cq_s[oilCompIdx];
298 const Value cq_sGas = cq_s[gasCompIdx];
299 const Value dis_gas = rs * cq_sOil;
300 const Value vap_oil = rv * cq_sGas;
302 cq_s[gasCompIdx] += dis_gas;
303 cq_s[oilCompIdx] += vap_oil;
306 if (this->isProducer()) {
307 perf_dis_gas_rate = getValue(dis_gas);
308 perf_vap_oil_rate = getValue(vap_oil);
311 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
312 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
313 const Value vap_wat = rvw * cq_sGas;
314 cq_s[waterCompIdx] += vap_wat;
315 if (this->isProducer())
316 perf_vap_wat_rate = getValue(vap_wat);
318 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
319 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
320 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
321 const Value cq_sWat = cq_s[waterCompIdx];
322 const Value cq_sGas = cq_s[gasCompIdx];
323 const Value vap_wat = rvw * cq_sGas;
324 const Value dis_gas_wat = rsw * cq_sWat;
325 cq_s[waterCompIdx] += vap_wat;
326 cq_s[gasCompIdx] += dis_gas_wat;
327 if (this->isProducer()) {
328 perf_vap_wat_rate = getValue(vap_wat);
329 perf_dis_gas_rate_in_water = getValue(dis_gas_wat);
335 if (!allow_cf && this->isProducer()) {
340 Value total_mob_dense = mob[0];
341 for (
int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
342 total_mob_dense += mob[componentIdx];
346 const Value cqt_i = - Tw * (total_mob_dense * drawdown);
349 Value volumeRatio = bhp * 0.0;
351 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
352 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
353 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
355 const Value d = 1.0 - rvw * rsw;
358 std::ostringstream sstr;
359 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
360 <<
" during computePerfRate calculations with rsw " << rsw
361 <<
", rvw " << rvw <<
" and pressure " << pressure
362 <<
" obtaining d " << d
363 <<
" Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) "
364 <<
" for this connection.";
365 deferred_logger.debug(sstr.str());
367 const Value tmp_wat = d > 0.0? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d : cmix_s[waterCompIdx];
368 volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
370 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d : cmix_s[waterCompIdx];
371 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
373 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
374 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
375 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
379 if constexpr (Indices::enableSolvent) {
380 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
383 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
384 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
385 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
387 const Value d = 1.0 - rv * rs;
390 std::ostringstream sstr;
391 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
392 <<
" during computePerfRate calculations with rs " << rs
393 <<
", rv " << rv <<
" and pressure " << pressure
394 <<
" obtaining d " << d
395 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
396 <<
" for this connection.";
397 deferred_logger.debug(sstr.str());
399 const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
400 volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
402 const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
403 volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
406 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
407 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
408 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
410 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
411 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
412 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
417 Value cqt_is = cqt_i/volumeRatio;
418 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
419 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
423 if (this->isProducer()) {
424 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
425 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
426 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
435 const double d = 1.0 - getValue(rv) * getValue(rs);
438 std::ostringstream sstr;
439 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
440 <<
" during computePerfRate calculations with rs " << rs
441 <<
", rv " << rv <<
" and pressure " << pressure
442 <<
" obtaining d " << d
443 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
444 <<
" for this connection.";
445 deferred_logger.debug(sstr.str());
449 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
452 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
454 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
459 perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
462 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
464 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
465 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
466 perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]);
467 perf_dis_gas_rate_in_water = getValue(rsw) * getValue(cq_s[waterCompIdx]);
474 template<
typename TypeTag>
476 StandardWell<TypeTag>::
477 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
479 const Well::InjectionControls& inj_controls,
480 const Well::ProductionControls& prod_controls,
481 WellState& well_state,
482 const GroupState& group_state,
483 DeferredLogger& deferred_logger)
487 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
490 this->linSys_.clear();
492 assembleWellEqWithoutIterationImpl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
498 template<
typename TypeTag>
500 StandardWell<TypeTag>::
501 assembleWellEqWithoutIterationImpl(
const Simulator& ebosSimulator,
503 const Well::InjectionControls& inj_controls,
504 const Well::ProductionControls& prod_controls,
505 WellState& well_state,
506 const GroupState& group_state,
507 DeferredLogger& deferred_logger)
510 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
511 const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
513 auto& ws = well_state.well(this->index_of_well_);
515 ws.vaporized_oil_rate = 0;
516 ws.dissolved_gas_rate = 0;
517 ws.dissolved_gas_rate_in_water = 0;
518 ws.vaporized_wat_rate = 0;
520 const int np = this->number_of_phases_;
522 std::vector<RateVector> connectionRates = this->connectionRates_;
523 auto& perf_data = ws.perf_data;
524 auto& perf_rates = perf_data.phase_rates;
525 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
527 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
528 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
529 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
530 calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
533 if constexpr (has_polymer && Base::has_polymermw) {
534 if (this->isInjector()) {
535 handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
538 const int cell_idx = this->well_cells_[perf];
539 for (
int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
541 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
543 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
545 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
546 assemblePerforationEq(cq_s_effective,
549 this->primary_variables_.numWellEq(),
553 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
554 auto& perf_rate_solvent = perf_data.solvent_rates;
555 perf_rate_solvent[perf] = cq_s[componentIdx].value();
557 perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
561 if constexpr (has_zFraction) {
562 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
563 assembleZFracEq(cq_s_zfrac_effective,
565 this->primary_variables_.numWellEq(),
570 this->connectionRates_ = connectionRates;
575 const auto& comm = this->parallel_well_info_.communication();
576 ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
577 ws.dissolved_gas_rate_in_water = comm.sum(ws.dissolved_gas_rate_in_water);
578 ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
579 ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
583 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
586 for (
int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
589 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
590 if (FluidSystem::numActivePhases() > 1) {
592 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
593 this->F0_[componentIdx]) * volume / dt;
595 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
596 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
597 assembleSourceEq(resWell_loc,
599 this->primary_variables_.numWellEq(),
603 const auto& summaryState = ebosSimulator.vanguard().summaryState();
604 const Schedule& schedule = ebosSimulator.vanguard().schedule();
605 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
606 assembleControlEq(well_state, group_state,
607 schedule, summaryState,
608 inj_controls, prod_controls,
609 this->primary_variables_,
610 this->connections_.rho(),
617 this->linSys_.invert();
619 OPM_DEFLOG_THROW(NumericalProblem,
"Error when inverting local well equations for well " + name(), deferred_logger);
626 template<
typename TypeTag>
628 StandardWell<TypeTag>::
629 calculateSinglePerf(
const Simulator& ebosSimulator,
631 WellState& well_state,
632 std::vector<RateVector>& connectionRates,
633 std::vector<EvalWell>& cq_s,
634 EvalWell& water_flux_s,
635 EvalWell& cq_s_zfrac_effective,
636 DeferredLogger& deferred_logger)
const
638 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
639 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
640 const int cell_idx = this->well_cells_[perf];
641 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
642 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
643 getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
645 double perf_dis_gas_rate = 0.;
646 double perf_dis_gas_rate_in_water = 0.;
647 double perf_vap_oil_rate = 0.;
648 double perf_vap_wat_rate = 0.;
649 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
650 const double Tw = this->well_index_[perf] * trans_mult;
651 computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
652 cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
654 auto& ws = well_state.well(this->index_of_well_);
655 auto& perf_data = ws.perf_data;
656 if constexpr (has_polymer && Base::has_polymermw) {
657 if (this->isInjector()) {
660 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
661 water_flux_s = cq_s[water_comp_idx];
664 handleInjectivityRate(ebosSimulator, perf, cq_s);
669 if (this->isProducer()) {
670 ws.dissolved_gas_rate += perf_dis_gas_rate;
671 ws.dissolved_gas_rate_in_water += perf_dis_gas_rate_in_water;
672 ws.vaporized_oil_rate += perf_vap_oil_rate;
673 ws.vaporized_wat_rate += perf_vap_wat_rate;
676 if constexpr (has_energy) {
677 connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
680 if constexpr (has_energy) {
682 auto fs = intQuants.fluidState();
683 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
684 if (!FluidSystem::phaseIsActive(phaseIdx)) {
689 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
690 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
691 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
692 if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
693 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
696 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
697 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
702 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
704 std::ostringstream sstr;
705 sstr <<
"Problematic d value " << d <<
" obtained for well " << this->name()
706 <<
" during calculateSinglePerf with rs " << fs.Rs()
707 <<
", rv " << fs.Rv()
708 <<
" obtaining d " << d
709 <<
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
710 <<
" for this connection.";
711 deferred_logger.debug(sstr.str());
712 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
714 if(FluidSystem::gasPhaseIdx == phaseIdx) {
715 cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
716 }
else if(FluidSystem::oilPhaseIdx == phaseIdx) {
718 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
724 if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
726 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
727 fs.setTemperature(this->well_ecl_.temperature());
728 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
729 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
730 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
731 paramCache.setRegionIndex(pvtRegionIdx);
732 paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
733 paramCache.updatePhase(fs, phaseIdx);
735 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
736 fs.setDensity(phaseIdx, rho);
737 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
738 fs.setEnthalpy(phaseIdx, h);
739 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
740 connectionRates[perf][Indices::contiEnergyEqIdx] += getValue(cq_r_thermal);
743 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
744 connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
749 if constexpr (has_polymer) {
751 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
752 EvalWell cq_s_poly = cq_s[waterCompIdx];
753 if (this->isInjector()) {
754 cq_s_poly *= this->wpolymer();
756 cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
759 auto& perf_rate_polymer = perf_data.polymer_rates;
760 perf_rate_polymer[perf] = cq_s_poly.value();
762 cq_s_poly *= this->well_efficiency_factor_;
763 connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
765 if constexpr (Base::has_polymermw) {
766 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
770 if constexpr (has_foam) {
772 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
773 EvalWell cq_s_foam = cq_s[gasCompIdx] * this->well_efficiency_factor_;
774 if (this->isInjector()) {
775 cq_s_foam *= this->wfoam();
777 cq_s_foam *= this->extendEval(intQuants.foamConcentration());
779 connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
782 if constexpr (has_zFraction) {
784 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
785 cq_s_zfrac_effective = cq_s[gasCompIdx];
786 if (this->isInjector()) {
787 cq_s_zfrac_effective *= this->wsolvent();
788 }
else if (cq_s_zfrac_effective.value() != 0.0) {
789 const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
790 cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
792 auto& perf_rate_solvent = perf_data.solvent_rates;
793 perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
795 cq_s_zfrac_effective *= this->well_efficiency_factor_;
796 connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
799 if constexpr (has_brine) {
801 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
803 EvalWell cq_s_sm = cq_s[waterCompIdx] - perf_vap_wat_rate;
804 if (this->isInjector()) {
805 cq_s_sm *= this->wsalt();
807 cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration());
810 auto& perf_rate_brine = perf_data.brine_rates;
811 perf_rate_brine[perf] = cq_s_sm.value();
813 cq_s_sm *= this->well_efficiency_factor_;
814 connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
817 if constexpr (has_micp) {
818 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
819 EvalWell cq_s_microbe = cq_s[waterCompIdx];
820 if (this->isInjector()) {
821 cq_s_microbe *= this->wmicrobes();
823 cq_s_microbe *= this->extendEval(intQuants.microbialConcentration());
825 connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe);
826 EvalWell cq_s_oxygen = cq_s[waterCompIdx];
827 if (this->isInjector()) {
828 cq_s_oxygen *= this->woxygen();
830 cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration());
832 connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen);
833 EvalWell cq_s_urea = cq_s[waterCompIdx];
834 if (this->isInjector()) {
835 cq_s_urea *= this->wurea();
837 cq_s_urea *= this->extendEval(intQuants.ureaConcentration());
839 connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea);
843 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
849 template<
typename TypeTag>
851 StandardWell<TypeTag>::
852 getMobilityEval(
const Simulator& ebosSimulator,
854 std::vector<EvalWell>& mob,
855 DeferredLogger& deferred_logger)
const
857 const int cell_idx = this->well_cells_[perf];
858 assert (
int(mob.size()) == this->num_components_);
859 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
860 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
864 const int satid = this->saturation_table_number_[perf] - 1;
865 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
866 if( satid == satid_elem ) {
868 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
869 if (!FluidSystem::phaseIsActive(phaseIdx)) {
873 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
874 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
877 mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
881 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
882 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
883 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
886 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
889 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
890 if (!FluidSystem::phaseIsActive(phaseIdx)) {
894 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
895 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
899 if constexpr (has_solvent) {
900 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
905 if constexpr (has_polymer) {
906 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
907 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
912 if constexpr (!Base::has_polymermw) {
913 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
918 template<
typename TypeTag>
920 StandardWell<TypeTag>::
921 getMobilityScalar(
const Simulator& ebosSimulator,
923 std::vector<Scalar>& mob,
924 DeferredLogger& deferred_logger)
const
926 const int cell_idx = this->well_cells_[perf];
927 assert (
int(mob.size()) == this->num_components_);
928 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
929 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
933 const int satid = this->saturation_table_number_[perf] - 1;
934 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
935 if( satid == satid_elem ) {
937 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
938 if (!FluidSystem::phaseIsActive(phaseIdx)) {
942 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
943 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
946 mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
950 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
951 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
952 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
955 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
958 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
959 if (!FluidSystem::phaseIsActive(phaseIdx)) {
963 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
964 mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx));
968 if constexpr (has_solvent) {
969 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
974 if constexpr (has_polymer) {
975 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
976 OPM_DEFLOG_THROW(std::runtime_error,
"Water is required when polymer is active", deferred_logger);
981 if constexpr (!Base::has_polymermw) {
982 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
983 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
984 for (
size_t i = 0; i < mob.size(); ++i) {
985 mob[i] = getValue(mob_eval[i]);
993 template<
typename TypeTag>
995 StandardWell<TypeTag>::
996 updateWellState(
const SummaryState& summary_state,
997 const BVectorWell& dwells,
998 WellState& well_state,
999 DeferredLogger& deferred_logger)
1001 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1003 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
1004 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
1006 updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
1007 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
1014 template<
typename TypeTag>
1016 StandardWell<TypeTag>::
1017 updatePrimaryVariablesNewton(
const BVectorWell& dwells,
1018 const bool stop_or_zero_rate_target,
1019 DeferredLogger& deferred_logger)
1021 const double dFLimit = this->param_.dwell_fraction_max_;
1022 const double dBHPLimit = this->param_.dbhp_max_rel_;
1023 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit);
1026 if constexpr (Base::has_polymermw) {
1027 this->primary_variables_.updateNewtonPolyMW(dwells);
1030 this->primary_variables_.checkFinite(deferred_logger);
1037 template<
typename TypeTag>
1039 StandardWell<TypeTag>::
1040 updateWellStateFromPrimaryVariables(
const bool stop_or_zero_rate_target,
1041 WellState& well_state,
1042 DeferredLogger& deferred_logger)
const
1044 this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
1047 if constexpr (Base::has_polymermw) {
1048 this->primary_variables_.copyToWellStatePolyMW(well_state);
1056 template<
typename TypeTag>
1058 StandardWell<TypeTag>::
1059 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1064 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1065 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1067 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1068 std::vector<Scalar> mob(this->num_components_, 0.0);
1069 getMobilityScalar(ebos_simulator, perf, mob, deferred_logger);
1071 const int cell_idx = this->well_cells_[perf];
1072 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1073 const auto& fs = int_quantities.fluidState();
1075 double p_r = this->getPerfCellPressure(fs).value();
1078 std::vector<double> b_perf(this->num_components_);
1079 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1080 if (!FluidSystem::phaseIsActive(phase)) {
1083 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1084 b_perf[comp_idx] = fs.invB(phase).value();
1086 if constexpr (has_solvent) {
1087 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
1091 const double h_perf = this->connections_.pressure_diff(perf);
1092 const double pressure_diff = p_r - h_perf;
1097 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1098 deferred_logger.debug(
"CROSSFLOW_IPR",
1099 "cross flow found when updateIPR for well " + name()
1100 +
" . The connection is ignored in IPR calculations");
1106 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1108 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1109 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1110 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1111 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1112 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1113 ipr_b_perf[comp_idx] += tw_mob;
1117 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1118 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1119 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1120 const double rs = (fs.Rs()).value();
1121 const double rv = (fs.Rv()).value();
1123 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1124 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1126 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1127 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1129 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1130 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1132 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1133 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1136 for (
size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1137 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1138 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1141 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1142 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1146 template<
typename TypeTag>
1148 StandardWell<TypeTag>::
1149 checkOperabilityUnderBHPLimit(
const WellState& well_state,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1151 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1152 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1155 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1156 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1159 double total_ipr_mass_rate = 0.0;
1160 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1162 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1166 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1167 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1169 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1170 total_ipr_mass_rate += ipr_rate * rho;
1172 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1173 this->operability_status_.operable_under_only_bhp_limit =
false;
1177 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1181 std::vector<double> well_rates_bhp_limit;
1182 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1184 this->adaptRatesForVFP(well_rates_bhp_limit);
1185 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1187 this->connections_.rho(),
1188 this->getALQ(well_state),
1190 const double thp_limit = this->getTHPConstraint(summaryState);
1191 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1192 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1203 this->operability_status_.operable_under_only_bhp_limit =
true;
1204 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1212 template<
typename TypeTag>
1214 StandardWell<TypeTag>::
1215 checkOperabilityUnderTHPLimit(
const Simulator& ebos_simulator,
const WellState& well_state, DeferredLogger& deferred_logger)
1217 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1218 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger)
1219 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1222 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1224 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1225 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1226 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1228 const double thp_limit = this->getTHPConstraint(summaryState);
1229 if (this->isProducer() && *obtain_bhp < thp_limit) {
1230 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1231 +
" bars is SMALLER than thp limit "
1232 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1233 +
" bars as a producer for well " + name();
1234 deferred_logger.debug(msg);
1236 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1237 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1238 +
" bars is LARGER than thp limit "
1239 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1240 +
" bars as a injector for well " + name();
1241 deferred_logger.debug(msg);
1244 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1245 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1246 if (!this->wellIsStopped()) {
1247 const double thp_limit = this->getTHPConstraint(summaryState);
1248 deferred_logger.debug(
" could not find bhp value at thp limit "
1249 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1250 +
" bar for well " + name() +
", the well might need to be closed ");
1259 template<
typename TypeTag>
1261 StandardWell<TypeTag>::
1262 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1264 bool all_drawdown_wrong_direction =
true;
1266 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1267 const int cell_idx = this->well_cells_[perf];
1268 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1269 const auto& fs = intQuants.fluidState();
1271 const double pressure = this->getPerfCellPressure(fs).value();
1272 const double bhp = this->primary_variables_.eval(Bhp).value();
1275 const double well_pressure = bhp + this->connections_.pressure_diff(perf);
1276 const double drawdown = pressure - well_pressure;
1281 if ( (drawdown < 0. && this->isInjector()) ||
1282 (drawdown > 0. && this->isProducer()) ) {
1283 all_drawdown_wrong_direction =
false;
1288 const auto& comm = this->parallel_well_info_.communication();
1289 if (comm.size() > 1)
1291 all_drawdown_wrong_direction =
1292 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1295 return all_drawdown_wrong_direction;
1301 template<
typename TypeTag>
1303 StandardWell<TypeTag>::
1304 canProduceInjectWithCurrentBhp(
const Simulator& ebos_simulator,
1305 const WellState& well_state,
1306 DeferredLogger& deferred_logger)
1308 const double bhp = well_state.well(this->index_of_well_).bhp;
1309 std::vector<double> well_rates;
1310 computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
1312 const double sign = (this->isProducer()) ? -1. : 1.;
1313 const double threshold = sign * std::numeric_limits<double>::min();
1315 bool can_produce_inject =
false;
1316 for (
const auto value : well_rates) {
1317 if (this->isProducer() && value < threshold) {
1318 can_produce_inject =
true;
1320 }
else if (this->isInjector() && value > threshold) {
1321 can_produce_inject =
true;
1326 if (!can_produce_inject) {
1327 deferred_logger.debug(
" well " + name() +
" CANNOT produce or inejct ");
1330 return can_produce_inject;
1337 template<
typename TypeTag>
1339 StandardWell<TypeTag>::
1340 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1342 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1348 template<
typename TypeTag>
1350 StandardWell<TypeTag>::
1351 computePropertiesForWellConnectionPressures(
const Simulator& ebosSimulator,
1352 const WellState& well_state,
1353 std::vector<double>& b_perf,
1354 std::vector<double>& rsmax_perf,
1355 std::vector<double>& rvmax_perf,
1356 std::vector<double>& rvwmax_perf,
1357 std::vector<double>& rswmax_perf,
1358 std::vector<double>& surf_dens_perf)
const
1360 std::function<Scalar(
int,
int)> getTemperature =
1361 [&ebosSimulator](
int cell_idx,
int phase_idx)
1363 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().temperature(phase_idx).value();
1365 std::function<Scalar(
int)> getSaltConcentration =
1366 [&ebosSimulator](
int cell_idx)
1368 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().saltConcentration().value();
1370 std::function<int(
int)> getPvtRegionIdx =
1371 [&ebosSimulator](
int cell_idx)
1373 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().pvtRegionIndex();
1375 std::function<Scalar(
int)> getInvFac =
1376 [&ebosSimulator](
int cell_idx)
1378 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1380 std::function<Scalar(
int)> getSolventDensity =
1381 [&ebosSimulator](
int cell_idx)
1383 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventRefDensity();
1386 this->connections_.computePropertiesForPressures(well_state,
1388 getSaltConcentration,
1404 template<
typename TypeTag>
1408 const std::vector<double>&
B_avg,
1414 assert((
int(
B_avg.size()) ==
this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1416 std::vector<double>
res;
1419 this->param_.max_residual_allowed_,
1420 this->param_.tolerance_wells_,
1421 this->param_.relaxed_tolerance_flow_well_,
1425 checkConvergenceExtraEqs(
res, report);
1434 template<
typename TypeTag>
1445 return ebosSimulator.model()
1446 .intensiveQuantities(
cell_idx, 0).fluidState();
1449 const int np = this->number_of_phases_;
1450 auto setToZero = [np](
double* x) ->
void
1452 std::fill_n(x, np, 0.0);
1455 auto addVector = [np](
const double* src,
double* dest) ->
void
1457 std::transform(src, src + np, dest, dest, std::plus<>{});
1460 auto& ws = well_state.well(this->index_of_well_);
1461 auto& perf_data = ws.perf_data;
1462 auto* wellPI = ws.productivity_index.data();
1463 auto* connPI = perf_data.prod_index.data();
1467 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1468 auto subsetPerfID = 0;
1470 for (
const auto& perf : *this->perf_data_) {
1471 auto allPerfID = perf.ecl_index;
1473 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
1478 std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
1479 getMobilityEval(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
1481 const auto& fs = fluidState(subsetPerfID);
1484 if (this->isInjector()) {
1485 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1486 mob, connPI, deferred_logger);
1489 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1492 addVector(connPI, wellPI);
1499 const auto& comm = this->parallel_well_info_.communication();
1500 if (comm.size() > 1) {
1501 comm.sum(wellPI, np);
1504 assert ((
static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1505 "Internal logic error in processing connections for PI/II");
1510 template<
typename TypeTag>
1512 StandardWell<TypeTag>::
1513 computeWellConnectionDensitesPressures(
const Simulator& ebosSimulator,
1514 const WellState& well_state,
1515 const std::vector<double>& b_perf,
1516 const std::vector<double>& rsmax_perf,
1517 const std::vector<double>& rvmax_perf,
1518 const std::vector<double>& rvwmax_perf,
1519 const std::vector<double>& rswmax_perf,
1520 const std::vector<double>& surf_dens_perf,
1521 DeferredLogger& deferred_logger)
1523 std::function<Scalar(
int,
int)> invB =
1524 [&ebosSimulator](
int cell_idx,
int phase_idx)
1526 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).fluidState().invB(phase_idx).value();
1528 std::function<Scalar(
int,
int)> mobility =
1529 [&ebosSimulator](
int cell_idx,
int phase_idx)
1531 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).mobility(phase_idx).value();
1533 std::function<Scalar(
int)> invFac =
1534 [&ebosSimulator](
int cell_idx)
1536 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
1538 std::function<Scalar(
int)> solventMobility =
1539 [&ebosSimulator](
int cell_idx)
1541 return ebosSimulator.model().intensiveQuantities(cell_idx, 0).solventMobility().value();
1544 this->connections_.computeProperties(well_state,
1562 template<
typename TypeTag>
1564 StandardWell<TypeTag>::
1565 computeWellConnectionPressures(
const Simulator& ebosSimulator,
1566 const WellState& well_state,
1567 DeferredLogger& deferred_logger)
1572 std::vector<double> b_perf;
1573 std::vector<double> rsmax_perf;
1574 std::vector<double> rvmax_perf;
1575 std::vector<double> rvwmax_perf;
1576 std::vector<double> rswmax_perf;
1577 std::vector<double> surf_dens_perf;
1578 computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf);
1579 computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf, deferred_logger);
1586 template<
typename TypeTag>
1588 StandardWell<TypeTag>::
1589 solveEqAndUpdateWellState(
const SummaryState& summary_state,
1590 WellState& well_state,
1591 DeferredLogger& deferred_logger)
1593 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1597 BVectorWell dx_well(1);
1598 dx_well[0].resize(this->primary_variables_.numWellEq());
1599 this->linSys_.solve( dx_well);
1601 updateWellState(summary_state, dx_well, well_state, deferred_logger);
1608 template<
typename TypeTag>
1610 StandardWell<TypeTag>::
1611 calculateExplicitQuantities(
const Simulator& ebosSimulator,
1612 const WellState& well_state,
1613 DeferredLogger& deferred_logger)
1615 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1616 updatePrimaryVariables(summary_state, well_state, deferred_logger);
1617 initPrimaryVariablesEvaluation();
1618 computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
1619 this->computeAccumWell();
1624 template<
typename TypeTag>
1627 apply(
const BVector& x, BVector&
Ax)
const
1629 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1631 if (this->param_.matrix_add_well_contributions_)
1643 template<
typename TypeTag>
1648 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1656 template<
typename TypeTag>
1664 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1667 xw[0].resize(this->primary_variables_.numWellEq());
1669 this->linSys_.recoverSolutionWell(x,
xw);
1676 template<
typename TypeTag>
1685 const int np = this->number_of_phases_;
1688 const bool allow_cf = this->getAllowCrossFlow();
1690 for (
int perf = 0;
perf < this->number_of_perforations_; ++
perf) {
1694 std::vector<Scalar>
mob(this->num_components_, 0.);
1699 std::vector<Scalar>
cq_s(this->num_components_, 0.);
1703 for(
int p = 0;
p <
np; ++
p) {
1707 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1712 template<
typename TypeTag>
1714 StandardWell<TypeTag>::
1715 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
1717 std::vector<double>& well_flux,
1718 DeferredLogger& deferred_logger)
const
1722 StandardWell<TypeTag> well_copy(*
this);
1727 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1728 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1731 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1732 auto inj_controls = well_copy.well_ecl_.isInjector()
1733 ? well_copy.well_ecl_.injectionControls(summary_state)
1734 : Well::InjectionControls(0);
1735 auto prod_controls = well_copy.well_ecl_.isProducer()
1736 ? well_copy.well_ecl_.productionControls(summary_state) :
1737 Well::ProductionControls(0);
1740 auto& ws = well_state_copy.well(this->index_of_well_);
1741 if (well_copy.well_ecl_.isInjector()) {
1742 inj_controls.bhp_limit = bhp;
1743 ws.injection_cmode = Well::InjectorCMode::BHP;
1745 prod_controls.bhp_limit = bhp;
1746 ws.production_cmode = Well::ProducerCMode::BHP;
1751 const int np = this->number_of_phases_;
1752 const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1753 for (
int phase = 0; phase < np; ++phase){
1754 well_state_copy.wellRates(this->index_of_well_)[phase]
1755 = sign * ws.well_potentials[phase];
1757 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1759 const double dt = ebosSimulator.timeStepSize();
1760 const bool converged = well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1762 const std::string msg =
" well " + name() +
" did not get converged during well potential calculations "
1763 " potentials are computed based on unconverged solution";
1764 deferred_logger.debug(msg);
1766 well_copy.updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
1767 well_copy.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
1768 well_copy.initPrimaryVariablesEvaluation();
1769 well_copy.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1775 template<
typename TypeTag>
1777 StandardWell<TypeTag>::
1778 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
1779 DeferredLogger& deferred_logger,
1780 const WellState &well_state)
const
1782 std::vector<double> potentials(this->number_of_phases_, 0.0);
1783 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1785 const auto& well = this->well_ecl_;
1786 if (well.isInjector()){
1787 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1788 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1789 if (bhp_at_thp_limit) {
1790 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1791 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1793 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1794 "Failed in getting converged thp based potential calculation for well "
1795 + name() +
". Instead the bhp based value is used");
1796 const double bhp = controls.bhp_limit;
1797 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1800 computeWellRatesWithThpAlqProd(
1801 ebos_simulator, summary_state,
1802 deferred_logger, potentials, this->getALQ(well_state)
1811 template<
typename TypeTag>
1813 StandardWell<TypeTag>::
1814 updateWellStateWithTHPTargetProd(
const Simulator& ebos_simulator,
1815 WellState& well_state,
1816 DeferredLogger& deferred_logger)
const
1818 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1820 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1821 ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
1822 if (bhp_at_thp_limit) {
1823 std::vector<double> rates(this->number_of_phases_, 0.0);
1824 computeWellRatesWithBhp(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
1825 auto& ws = well_state.well(this->name());
1826 ws.surface_rates = rates;
1827 ws.bhp = *bhp_at_thp_limit;
1828 ws.thp = this->getTHPConstraint(summary_state);
1837 template<
typename TypeTag>
1839 StandardWell<TypeTag>::
1840 computeWellRatesAndBhpWithThpAlqProd(
const Simulator &ebos_simulator,
1841 const SummaryState &summary_state,
1842 DeferredLogger &deferred_logger,
1843 std::vector<double> &potentials,
1847 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1848 ebos_simulator, summary_state, alq, deferred_logger);
1849 if (bhp_at_thp_limit) {
1850 const auto& controls = this->well_ecl_.productionControls(summary_state);
1851 bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1852 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1855 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
1856 "Failed in getting converged thp based potential calculation for well "
1857 + name() +
". Instead the bhp based value is used");
1858 const auto& controls = this->well_ecl_.productionControls(summary_state);
1859 bhp = controls.bhp_limit;
1860 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1865 template<
typename TypeTag>
1867 StandardWell<TypeTag>::
1868 computeWellRatesWithThpAlqProd(
const Simulator &ebos_simulator,
1869 const SummaryState &summary_state,
1870 DeferredLogger &deferred_logger,
1871 std::vector<double> &potentials,
1875 computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1882 template<
typename TypeTag>
1887 std::vector<double>& well_potentials,
1890 const int np = this->number_of_phases_;
1891 well_potentials.resize(
np, 0.0);
1893 if (this->wellIsStopped()) {
1897 this->operability_status_.has_negative_potentials =
false;
1901 const auto&
ws = well_state.well(this->index_of_well_);
1902 if (this->isInjector()) {
1903 const Well::InjectorCMode&
current =
ws.injection_cmode;
1904 if (
current == Well::InjectorCMode::THP) {
1907 if (
current == Well::InjectorCMode::BHP) {
1911 const Well::ProducerCMode&
current =
ws.production_cmode;
1912 if (
current == Well::ProducerCMode::THP) {
1915 if (
current == Well::ProducerCMode::BHP) {
1922 const double sign = this->isInjector() ? 1.0:-1.0;
1923 for (
int phase = 0; phase <
np; ++phase){
1930 for (
int phase = 0; phase <
np; ++phase){
1931 well_potentials[phase] =
sign *
ws.surface_rates[phase];
1938 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1949 if (this->isInjector())
1950 bhp = std::max(
ws.bhp, bhp);
1952 bhp = std::min(
ws.bhp, bhp);
1954 assert(std::abs(bhp) != std::numeric_limits<double>::max());
1955 computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials,
deferred_logger);
1958 well_potentials = computeWellPotentialWithTHP(ebosSimulator,
deferred_logger, well_state);
1961 const double sign = this->isInjector() ? 1.0:-1.0;
1963 for (
int phase = 0; phase <
np; ++phase){
1964 well_potentials[phase] *=
sign;
1969 this->operability_status_.has_negative_potentials =
true;
1970 const std::string
msg = std::string(
"well ") + this->name() + std::string(
": has negative potentials and is not operable");
1979 template<
typename TypeTag>
1986 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1992 if constexpr (Base::has_polymermw) {
1993 this->primary_variables_.updatePolyMW(well_state);
1996 this->primary_variables_.checkFinite(deferred_logger);
2002 template<
typename TypeTag>
2004 StandardWell<TypeTag>::
2005 getRefDensity()
const
2007 return this->connections_.rho();
2013 template<
typename TypeTag>
2015 StandardWell<TypeTag>::
2016 updateWaterMobilityWithPolymer(
const Simulator& ebos_simulator,
2018 std::vector<EvalWell>& mob,
2019 DeferredLogger& deferred_logger)
const
2021 const int cell_idx = this->well_cells_[perf];
2022 const auto& int_quant = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
2023 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
2027 if (this->isInjector()) {
2029 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
2030 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2031 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration,
true) );
2034 if (PolymerModule::hasPlyshlog()) {
2037 if (this->isInjector() && this->wpolymer() == 0.) {
2042 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
2043 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2045 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
2046 double perf_dis_gas_rate = 0.;
2047 double perf_dis_gas_rate_in_water = 0.;
2048 double perf_vap_oil_rate = 0.;
2049 double perf_vap_wat_rate = 0.;
2050 double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
2051 const double Tw = this->well_index_[perf] * trans_mult;
2052 computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
2053 cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
2055 const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
2056 const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
2057 const auto& scaled_drainage_info =
2058 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
2059 const double swcr = scaled_drainage_info.Swcr;
2060 const EvalWell poro = this->extendEval(int_quant.porosity());
2061 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
2063 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
2064 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2065 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
2067 if (PolymerModule::hasShrate()) {
2070 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
2072 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
2073 int_quant.pvtRegionIndex(),
2076 mob[waterCompIdx] /= shear_factor;
2080 template<
typename TypeTag>
2082 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian)
const
2084 this->linSys_.extract(jacobian);
2088 template <
typename TypeTag>
2090 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
2091 const BVector& weights,
2092 const int pressureVarIndex,
2093 const bool use_well_weights,
2094 const WellState& well_state)
const
2096 this->linSys_.extractCPRPressureMatrix(jacobian,
2107 template<
typename TypeTag>
2108 typename StandardWell<TypeTag>::EvalWell
2109 StandardWell<TypeTag>::
2110 pskinwater(
const double throughput,
2111 const EvalWell& water_velocity,
2112 DeferredLogger& deferred_logger)
const
2114 if constexpr (Base::has_polymermw) {
2115 const int water_table_id = this->polymerWaterTable_();
2116 if (water_table_id <= 0) {
2117 OPM_DEFLOG_THROW(std::runtime_error,
2118 fmt::format(
"Unused SKPRWAT table id used for well {}", name()),
2121 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
2122 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2124 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2125 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
2128 OPM_DEFLOG_THROW(std::runtime_error,
2129 fmt::format(
"Polymermw is not activated, while injecting "
2130 "skin pressure is requested for well {}", name()),
2139 template<
typename TypeTag>
2140 typename StandardWell<TypeTag>::EvalWell
2141 StandardWell<TypeTag>::
2142 pskin(
const double throughput,
2143 const EvalWell& water_velocity,
2144 const EvalWell& poly_inj_conc,
2145 DeferredLogger& deferred_logger)
const
2147 if constexpr (Base::has_polymermw) {
2148 const double sign = water_velocity >= 0. ? 1.0 : -1.0;
2149 const EvalWell water_velocity_abs = abs(water_velocity);
2150 if (poly_inj_conc == 0.) {
2151 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
2153 const int polymer_table_id = this->polymerTable_();
2154 if (polymer_table_id <= 0) {
2155 OPM_DEFLOG_THROW(std::runtime_error,
2156 fmt::format(
"Unavailable SKPRPOLY table id used for well {}", name()),
2159 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
2160 const double reference_concentration = skprpolytable.refConcentration;
2161 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2163 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2164 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2165 if (poly_inj_conc == reference_concentration) {
2166 return sign * pskin_poly;
2169 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2170 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2171 return sign * pskin;
2173 OPM_DEFLOG_THROW(std::runtime_error,
2174 fmt::format(
"Polymermw is not activated, while injecting "
2175 "skin pressure is requested for well {}", name()),
2184 template<
typename TypeTag>
2185 typename StandardWell<TypeTag>::EvalWell
2186 StandardWell<TypeTag>::
2187 wpolymermw(
const double throughput,
2188 const EvalWell& water_velocity,
2189 DeferredLogger& deferred_logger)
const
2191 if constexpr (Base::has_polymermw) {
2192 const int table_id = this->polymerInjTable_();
2193 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2194 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2195 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2196 if (this->wpolymer() == 0.) {
2197 return molecular_weight;
2199 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2200 return molecular_weight;
2202 OPM_DEFLOG_THROW(std::runtime_error,
2203 fmt::format(
"Polymermw is not activated, while injecting "
2204 "polymer molecular weight is requested for well {}", name()),
2213 template<
typename TypeTag>
2215 StandardWell<TypeTag>::
2216 updateWaterThroughput(
const double dt, WellState &well_state)
const
2218 if constexpr (Base::has_polymermw) {
2219 if (this->isInjector()) {
2220 auto& ws = well_state.well(this->index_of_well_);
2221 auto& perf_water_throughput = ws.perf_data.water_throughput;
2222 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2223 const double perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
2225 if (perf_water_vel > 0.) {
2226 perf_water_throughput[perf] += perf_water_vel * dt;
2237 template<
typename TypeTag>
2239 StandardWell<TypeTag>::
2240 handleInjectivityRate(
const Simulator& ebosSimulator,
2242 std::vector<EvalWell>& cq_s)
const
2244 const int cell_idx = this->well_cells_[perf];
2245 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
2246 const auto& fs = int_quants.fluidState();
2247 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2248 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2249 const int wat_vel_index = Bhp + 1 + perf;
2250 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2254 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2260 template<
typename TypeTag>
2262 StandardWell<TypeTag>::
2263 handleInjectivityEquations(
const Simulator& ebosSimulator,
2264 const WellState& well_state,
2266 const EvalWell& water_flux_s,
2267 DeferredLogger& deferred_logger)
2269 const int cell_idx = this->well_cells_[perf];
2270 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
2271 const auto& fs = int_quants.fluidState();
2272 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2273 const EvalWell water_flux_r = water_flux_s / b_w;
2274 const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2275 const EvalWell water_velocity = water_flux_r / area;
2276 const int wat_vel_index = Bhp + 1 + perf;
2279 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2281 const auto& ws = well_state.well(this->index_of_well_);
2282 const auto& perf_data = ws.perf_data;
2283 const auto& perf_water_throughput = perf_data.water_throughput;
2284 const double throughput = perf_water_throughput[perf];
2285 const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2287 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2288 poly_conc.setValue(this->wpolymer());
2291 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2292 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2294 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
2295 assembleInjectivityEq(eq_pskin,
2300 this->primary_variables_.numWellEq(),
2308 template<
typename TypeTag>
2310 StandardWell<TypeTag>::
2311 checkConvergenceExtraEqs(
const std::vector<double>& res,
2312 ConvergenceReport& report)
const
2317 if constexpr (Base::has_polymermw) {
2318 WellConvergence(*this).
2319 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2327 template<
typename TypeTag>
2329 StandardWell<TypeTag>::
2330 updateConnectionRatePolyMW(
const EvalWell& cq_s_poly,
2331 const IntensiveQuantities& int_quants,
2332 const WellState& well_state,
2334 std::vector<RateVector>& connectionRates,
2335 DeferredLogger& deferred_logger)
const
2338 EvalWell cq_s_polymw = cq_s_poly;
2339 if (this->isInjector()) {
2340 const int wat_vel_index = Bhp + 1 + perf;
2341 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2342 if (water_velocity > 0.) {
2343 const auto& ws = well_state.well(this->index_of_well_);
2344 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2345 const double throughput = perf_water_throughput[perf];
2346 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2347 cq_s_polymw *= molecular_weight;
2353 }
else if (this->isProducer()) {
2354 if (cq_s_polymw < 0.) {
2355 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2362 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2370 template<
typename TypeTag>
2371 std::optional<double>
2372 StandardWell<TypeTag>::
2373 computeBhpAtThpLimitProd(
const WellState& well_state,
2374 const Simulator& ebos_simulator,
2375 const SummaryState& summary_state,
2376 DeferredLogger& deferred_logger)
const
2378 return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
2380 this->getALQ(well_state),
2384 template<
typename TypeTag>
2385 std::optional<double>
2386 StandardWell<TypeTag>::
2387 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
2388 const SummaryState& summary_state,
2389 const double alq_value,
2390 DeferredLogger& deferred_logger)
const
2393 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2399 std::vector<double> rates(3);
2400 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2401 this->adaptRatesForVFP(rates);
2405 double max_pressure = 0.0;
2406 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
2407 const int cell_idx = this->well_cells_[perf];
2408 const auto& int_quants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
2409 const auto& fs = int_quants.fluidState();
2410 double pressure_cell = this->getPerfCellPressure(fs).value();
2411 max_pressure = std::max(max_pressure, pressure_cell);
2413 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2416 this->connections_.rho(),
2418 this->getTHPConstraint(summary_state),
2422 auto v = frates(*bhpAtLimit);
2423 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2429 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2433 std::vector<double> rates(3);
2434 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
2435 this->adaptRatesForVFP(rates);
2439 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2442 this->connections_.rho(),
2444 this->getTHPConstraint(summary_state),
2450 auto v = frates(*bhpAtLimit);
2451 if (std::all_of(v.cbegin(), v.cend(), [](
double i){ return i <= 0; }) ) {
2457 return std::nullopt;
2462 template<
typename TypeTag>
2463 std::optional<double>
2464 StandardWell<TypeTag>::
2465 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
2466 const SummaryState& summary_state,
2467 DeferredLogger& deferred_logger)
const
2470 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
2476 std::vector<double> rates(3);
2477 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2481 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2483 this->connections_.rho(),
2494 template<
typename TypeTag>
2496 StandardWell<TypeTag>::
2497 iterateWellEqWithControl(
const Simulator& ebosSimulator,
2499 const Well::InjectionControls& inj_controls,
2500 const Well::ProductionControls& prod_controls,
2501 WellState& well_state,
2502 const GroupState& group_state,
2503 DeferredLogger& deferred_logger)
2505 const int max_iter = this->param_.max_inner_iter_wells_;
2508 bool relax_convergence =
false;
2509 this->regularize_ =
false;
2511 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2513 if (it > this->param_.strict_inner_iter_wells_) {
2514 relax_convergence =
true;
2515 this->regularize_ =
true;
2518 auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
2520 converged = report.converged();
2526 const auto& summary_state = ebosSimulator.vanguard().summaryState();
2527 solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
2534 initPrimaryVariablesEvaluation();
2535 }
while (it < max_iter);
2541 template<
typename TypeTag>
2548 std::vector<double>
well_q_s(this->num_components_, 0.);
2549 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2550 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2551 for (
int perf = 0;
perf < this->number_of_perforations_; ++
perf) {
2554 std::vector<Scalar>
mob(this->num_components_, 0.);
2556 std::vector<Scalar>
cq_s(this->num_components_, 0.);
2561 for (
int comp = 0;
comp < this->num_components_; ++
comp) {
2565 const auto& comm = this->parallel_well_info_.communication();
2566 if (comm.size() > 1)
2577 template <
typename TypeTag>
2581 const std::function<
double(
const double)>&
connPICalc,
2582 const std::vector<EvalWell>&
mobility,
2586 const int np = this->number_of_phases_;
2587 for (
int p = 0;
p <
np; ++
p) {
2591 mobility[ this->flowPhaseToEbosCompIdx(
p) ].value()
2592 *
fs.invB(this->flowPhaseToEbosPhaseIdx(
p)).value();
2597 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2598 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2600 const auto io = pu.phase_pos[Oil];
2601 const auto ig = pu.phase_pos[Gas];
2603 const auto vapoil = connPI[ig] * fs.Rv().value();
2604 const auto disgas = connPI[io] * fs.Rs().value();
2606 connPI[io] += vapoil;
2607 connPI[ig] += disgas;
2615 template <
typename TypeTag>
2617 StandardWell<TypeTag>::
2618 computeConnLevelInjInd(
const typename StandardWell<TypeTag>::FluidState& fs,
2619 const Phase preferred_phase,
2620 const std::function<
double(
const double)>& connIICalc,
2621 const std::vector<EvalWell>& mobility,
2623 DeferredLogger& deferred_logger)
const
2629 if (preferred_phase == Phase::GAS) {
2630 phase_pos = pu.phase_pos[Gas];
2632 else if (preferred_phase == Phase::OIL) {
2633 phase_pos = pu.phase_pos[Oil];
2635 else if (preferred_phase == Phase::WATER) {
2636 phase_pos = pu.phase_pos[Water];
2639 OPM_DEFLOG_THROW(NotImplemented,
2640 fmt::format(
"Unsupported Injector Type ({}) "
2641 "for well {} during connection I.I. calculation",
2642 static_cast<int>(preferred_phase), this->name()),
2646 const auto zero = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
2647 const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
2648 connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
Definition AquiferInterface.hpp:35
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
Definition StandardWell.hpp:60
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition StandardWell_impl.hpp:1627
Class for computing BHP limits.
Definition WellBhpThpCalculator.hpp:42
double mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
Definition WellBhpThpCalculator.cpp:85
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:60
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