22#ifndef OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
23#define OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
26#ifndef OPM_WELLINTERFACE_HEADER_INCLUDED
31#include <opm/common/Exceptions.hpp>
33#include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
34#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
43#include <dune/common/version.hh>
50#include <fmt/format.h>
56 template<
typename TypeTag>
63 const int pvtRegionIdx,
64 const int num_components,
66 const int index_of_well,
82 if (well.isInjector()) {
83 auto injectorType = this->
well_ecl_.injectorType();
84 if (injectorType == InjectorType::GAS) {
92 template<
typename TypeTag>
96 const std::vector<Scalar>& ,
98 const std::vector<Scalar>& B_avg,
99 const bool changed_to_open_this_step)
101 this->phase_usage_ = phase_usage_arg;
102 this->gravity_ = gravity_arg;
104 this->changed_to_open_this_step_ = changed_to_open_this_step;
110 template<
typename TypeTag>
115 if constexpr (has_polymer) {
116 return this->wpolymer_();
126 template<
typename TypeTag>
131 if constexpr (has_foam) {
132 return this->wfoam_();
140 template<
typename TypeTag>
145 if constexpr (has_brine) {
146 return this->wsalt_();
152 template<
typename TypeTag>
157 if constexpr (has_micp) {
158 return this->wmicrobes_();
164 template<
typename TypeTag>
169 if constexpr (has_micp) {
170 return this->woxygen_();
176 template<
typename TypeTag>
181 if constexpr (has_micp) {
182 return this->wurea_();
188 template<
typename TypeTag>
198 if (stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
202 const auto& summaryState = simulator.vanguard().summaryState();
203 const auto& schedule = simulator.vanguard().schedule();
204 const auto& well = this->well_ecl_;
205 auto& ws = well_state.
well(this->index_of_well_);
207 if (well.isInjector()) {
208 from = WellInjectorCMode2String(ws.injection_cmode);
210 from = WellProducerCMode2String(ws.production_cmode);
213 const int episodeIdx = simulator.episodeIndex();
214 const int iterationIdx = simulator.model().newtonMethod().numIterations();
215 const int nupcol = schedule[episodeIdx].nupcol();
216 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
219 const bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
221 const auto msg = fmt::format(
" The control mode for well {} is oscillating. \n"
222 "We don't allow for more than {} switches after NUPCOL iterations. (NUPCOL = {}) \n"
223 "The control is kept at {}.",
224 this->name(), this->param_.max_number_of_well_switches_, nupcol, from);
225 deferred_logger.
info(msg);
227 this->well_control_log_.push_back(from);
231 bool changed =
false;
232 if (iog == IndividualOrGroup::Individual) {
233 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
234 }
else if (iog == IndividualOrGroup::Group) {
235 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState,
true, deferred_logger);
237 assert(iog == IndividualOrGroup::Both);
238 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
244 if (well.isInjector()) {
245 to = WellInjectorCMode2String(ws.injection_cmode);
247 to = WellProducerCMode2String(ws.production_cmode);
249 std::ostringstream ss;
250 ss <<
" Switching control mode for well " << this->name()
254 ss <<
" on rank " << cc.rank();
256 deferred_logger.
debug(ss.str());
261 if (iterationIdx >= nupcol || this->well_control_log_.empty()) {
262 this->well_control_log_.push_back(from);
264 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
265 updatePrimaryVariables(simulator, well_state, deferred_logger);
271 template<
typename TypeTag>
277 const Well::InjectionControls& inj_controls,
278 const Well::ProductionControls& prod_controls,
281 const bool fixed_control,
282 const bool fixed_status)
285 const auto& summary_state = simulator.vanguard().summaryState();
286 const auto& schedule = simulator.vanguard().schedule();
287 auto& ws = well_state.
well(this->index_of_well_);
289 if (this->isInjector()) {
290 from = WellInjectorCMode2String(ws.injection_cmode);
292 from = WellProducerCMode2String(ws.production_cmode);
294 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
295 if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
299 const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
300 if (!this->wellIsStopped()){
301 if (wqTotal*sgn <= 0.0 && !fixed_status){
305 bool changed =
false;
306 if (!fixed_control) {
311 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
312 prod_controls.hasControl(Well::ProducerCMode::GRUP);
313 bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
314 if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
315 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
317 if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
318 changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,
false, deferred_logger);
322 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
323 ws.production_cmode == Well::ProducerCMode::THP;
325 ws.thp = this->getTHPConstraint(summary_state);
328 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
330 updatePrimaryVariables(simulator, well_state, deferred_logger);
335 }
else if (!fixed_status){
337 const Scalar bhp = well_state.
well(this->index_of_well_).bhp;
338 Scalar prod_limit = prod_controls.bhp_limit;
339 Scalar inj_limit = inj_controls.bhp_limit;
340 const bool has_thp = this->wellHasTHPConstraints(summary_state);
342 std::vector<Scalar> rates(this->num_components_);
343 if (this->isInjector()){
345 calculateBhpFromThp(well_state, rates,
348 this->getRefDensity(),
350 inj_limit = std::min(bhp_thp,
static_cast<Scalar>(inj_controls.bhp_limit));
355 calculateMinimumBhpFromThp(well_state,
358 this->getRefDensity());
359 prod_limit = std::max(bhp_min,
static_cast<Scalar>(prod_controls.bhp_limit));
362 const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
365 well_state.
well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
367 well_state.
well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
378 template<
typename TypeTag>
382 const double simulation_time,
385 WellTestState& well_test_state,
388 std::map<std::string, double>& open_times,
392 deferred_logger.
info(
" well " + this->name() +
" is being tested");
395 auto& ws = well_state_copy.
well(this->indexOfWell());
397 const auto& summary_state = simulator.vanguard().summaryState();
398 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
399 if (this->isProducer()) {
400 ws.production_cmode = has_thp_limit ? Well::ProducerCMode::THP : Well::ProducerCMode::BHP;
402 ws.injection_cmode = has_thp_limit ? Well::InjectorCMode::THP : Well::InjectorCMode::BHP;
407 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
408 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
409 updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
411 if (this->isProducer()) {
412 const auto& schedule = simulator.vanguard().schedule();
413 const auto report_step = simulator.episodeIndex();
414 const auto& glo = schedule.glo(report_step);
416 gliftBeginTimeStepWellTestUpdateALQ(simulator,
425 WellTestState welltest_state_temp;
427 bool testWell =
true;
432 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
433 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
435 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
436 deferred_logger.
debug(msg);
441 updateWellOperability(simulator, well_state_copy, deferred_logger);
442 if ( !this->isOperableAndSolvable() ) {
443 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
444 deferred_logger.
debug(msg);
447 std::vector<Scalar> potentials;
449 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
450 }
catch (
const std::exception& e) {
451 const std::string msg = fmt::format(
"well {}: computeWellPotentials() "
452 "failed during testing for re-opening: ",
453 this->name(), e.what());
454 deferred_logger.
info(msg);
457 const int np = well_state_copy.
numPhases();
458 for (
int p = 0; p < np; ++p) {
459 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
461 const bool under_zero_target = this->wellUnderZeroGroupRateTarget(simulator, well_state_copy, deferred_logger);
462 this->updateWellTestState(well_state_copy.
well(this->indexOfWell()),
468 this->closeCompletions(welltest_state_temp);
474 if ( welltest_state_temp.num_closed_wells() > 0 ||
475 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
481 if (!welltest_state_temp.well_is_closed(this->name())) {
482 well_test_state.open_well(this->name());
484 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
485 deferred_logger.
info(msg);
488 for (
const auto& completion : this->well_ecl_.getCompletions()) {
489 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
490 well_test_state.open_completion(this->name(), completion.first);
492 well_state = well_state_copy;
493 open_times.try_emplace(this->name(), well_test_state.lastTestTime(this->name()));
500 template<
typename TypeTag>
510 const auto& summary_state = simulator.vanguard().summaryState();
511 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
512 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
513 const auto& ws = well_state.
well(this->indexOfWell());
514 const auto pmode_orig = ws.production_cmode;
515 const auto imode_orig = ws.injection_cmode;
516 bool converged =
false;
519 if (!this->param_.local_well_solver_control_switching_){
520 converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
522 if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (this->well_ecl_.getStatus() == WellStatus::OPEN)) {
523 converged = solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
525 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
529 }
catch (NumericalProblem& e ) {
530 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
531 deferred_logger.
warning(
"INNER_ITERATION_FAILED", msg);
536 if (ws.production_cmode != pmode_orig || ws.injection_cmode != imode_orig) {
538 if (this->isInjector()) {
539 from = WellInjectorCMode2String(imode_orig);
540 to = WellInjectorCMode2String(ws.injection_cmode);
542 from = WellProducerCMode2String(pmode_orig);
543 to = WellProducerCMode2String(ws.production_cmode);
545 const auto msg = fmt::format(
" Well {} switched from {} to {} during local solve", this->name(), from, to);
546 deferred_logger.
debug(msg);
547 const int episodeIdx = simulator.episodeIndex();
548 const int iterationIdx = simulator.model().newtonMethod().numIterations();
549 const auto& schedule = simulator.vanguard().schedule();
550 const int nupcol = schedule[episodeIdx].nupcol();
554 if (iterationIdx >= nupcol || this->well_control_log_.empty()) {
555 this->well_control_log_.push_back(from);
563 template<
typename TypeTag>
568 const Well::InjectionControls& inj_controls,
569 const Well::ProductionControls& prod_controls,
575 const auto& summary_state = simulator.vanguard().summaryState();
576 bool is_operable =
true;
577 bool converged =
true;
578 auto& ws = well_state.
well(this->index_of_well_);
580 if (this->wellIsStopped()) {
582 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
583 if (!bhp_target.has_value()) {
585 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
586 deferred_logger.
debug(msg);
589 solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
593 ws.thp = this->getTHPConstraint(summary_state);
594 const Scalar bhp = std::max(bhp_target.value(),
595 static_cast<Scalar>(prod_controls.bhp_limit));
596 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
601 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
604 const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
606 if (converged && !stoppedOrZeroRateTarget(simulator, well_state, deferred_logger) && isThp) {
607 auto rates = well_state.
well(this->index_of_well_).surface_rates;
608 this->adaptRatesForVFP(rates);
609 this->updateIPRImplicit(simulator, well_state, deferred_logger);
613 this->operability_status_.use_vfpexplicit =
true;
616 const Scalar reltol = 1e-3;
617 const Scalar cur_bhp = ws.bhp;
618 if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
619 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
620 deferred_logger.
debug(msg);
621 solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger);
623 ws.thp = this->getTHPConstraint(summary_state);
624 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
631 this->operability_status_.use_vfpexplicit =
true;
633 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
634 if (!bhp_target.has_value()) {
638 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
642 const Scalar bhp = std::max(bhp_target.value(),
643 static_cast<Scalar>(prod_controls.bhp_limit));
644 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
645 ws.thp = this->getTHPConstraint(summary_state);
646 const auto msg = fmt::format(
"Well {} did not converge, re-solving with explicit fractions for VFP caculations.", this->name());
647 deferred_logger.
debug(msg);
648 converged = this->iterateWellEqWithSwitching(simulator, dt,
657 is_operable = is_operable && !this->wellIsStopped();
658 this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
659 this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
663 template<
typename TypeTag>
664 std::optional<typename WellInterface<TypeTag>::Scalar>
669 const SummaryState& summary_state,
672 if (!this->wellHasTHPConstraints(summary_state)) {
674 const bool converged = solveWellWithBhp(simulator, dt, bhp_limit, well_state, deferred_logger);
675 if (!converged || this->wellIsStopped()) {
686 const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
687 if (!converged || this->wellIsStopped()) {
690 this->updateIPRImplicit(simulator, well_state, deferred_logger);
691 auto rates = well_state.
well(this->index_of_well_).surface_rates;
692 this->adaptRatesForVFP(rates);
696 template<
typename TypeTag>
708 auto inj_controls = Well::InjectionControls(0);
709 auto prod_controls = Well::ProductionControls(0);
710 auto& ws = well_state.
well(this->index_of_well_);
711 auto cmode_inj = ws.injection_cmode;
712 auto cmode_prod = ws.production_cmode;
713 if (this->isInjector()) {
714 inj_controls.addControl(Well::InjectorCMode::BHP);
715 inj_controls.bhp_limit = bhp;
716 inj_controls.cmode = Well::InjectorCMode::BHP;
717 ws.injection_cmode = Well::InjectorCMode::BHP;
719 prod_controls.addControl(Well::ProducerCMode::BHP);
720 prod_controls.bhp_limit = bhp;
721 prod_controls.cmode = Well::ProducerCMode::BHP;
722 ws.production_cmode = Well::ProducerCMode::BHP;
727 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true);
728 ws.injection_cmode = cmode_inj;
729 ws.production_cmode = cmode_prod;
733 template<
typename TypeTag>
743 const auto well_status_orig = this->wellStatus_;
747 auto inj_controls = Well::InjectionControls(0);
748 auto prod_controls = Well::ProductionControls(0);
749 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true,
true);
750 this->wellStatus_ = well_status_orig;
754 template<
typename TypeTag>
763 const double dt = simulator.timeStepSize();
764 const bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
766 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" converged");
769 const int max_iter = this->param_.max_welleq_iter_;
770 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
776 template<
typename TypeTag>
785 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
790 const double dt = simulator.timeStepSize();
791 bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
801 auto& ws = well_state.
well(this->indexOfWell());
802 bool thp_control =
false;
803 if (this->well_ecl_.isInjector()) {
804 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
806 ws.injection_cmode = Well::InjectorCMode::BHP;
807 if (this->well_control_log_.empty()) {
808 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
812 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
814 ws.production_cmode = Well::ProducerCMode::BHP;
815 if (this->well_control_log_.empty()) {
816 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
821 const std::string msg = std::string(
"The newly opened well ") + this->name()
822 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
823 deferred_logger.
debug(msg);
824 converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
829 const int max_iter = this->param_.max_welleq_iter_;
830 deferred_logger.
debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
832 well_state = well_state0;
838 template <
typename TypeTag>
848 prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
849 assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
854 template <
typename TypeTag>
864 const auto& summary_state = simulator.vanguard().summaryState();
865 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
866 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
869 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
874 template<
typename TypeTag>
884 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
886 if (this->param_.check_well_operability_iter_)
887 checkWellOperability(simulator, well_state, deferred_logger);
890 const int iteration_idx = simulator.model().newtonMethod().numIterations();
891 if (iteration_idx < this->param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
892 const auto& ws = well_state.
well(this->indexOfWell());
893 const bool nonzero_rate_original =
894 std::any_of(ws.surface_rates.begin(),
895 ws.surface_rates.begin() + well_state.
numPhases(),
896 [](
Scalar rate) { return rate != Scalar(0.0); });
898 this->operability_status_.solvable =
true;
899 bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
902 const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
903 if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
907 this->operability_status_.resetOperability();
909 deferred_logger.
debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
913 if (this->param_.shut_unsolvable_wells_) {
914 this->operability_status_.solvable =
false;
918 if (this->operability_status_.has_negative_potentials) {
919 auto well_state_copy = well_state;
920 std::vector<Scalar> potentials;
922 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
923 }
catch (
const std::exception& e) {
924 const std::string msg = fmt::format(
"well {}: computeWellPotentials() failed "
925 "during attempt to recompute potentials for well: ",
926 this->name(), e.what());
927 deferred_logger.
info(msg);
928 this->operability_status_.has_negative_potentials =
true;
930 auto& ws = well_state.
well(this->indexOfWell());
932 for (
int p = 0; p < np; ++p) {
933 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
936 this->changed_to_open_this_step_ =
false;
937 const bool well_operable = this->operability_status_.isOperableAndSolvable();
939 if (!well_operable && old_well_operable) {
940 deferred_logger.
debug(
" well " + this->name() +
" gets STOPPED during iteration ");
942 changed_to_stopped_this_step_ =
true;
943 }
else if (well_operable && !old_well_operable) {
944 deferred_logger.
debug(
" well " + this->name() +
" gets REVIVED during iteration ");
946 changed_to_stopped_this_step_ =
false;
947 this->changed_to_open_this_step_ =
true;
951 template<
typename TypeTag>
955 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
958 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
959 if (this->cells()[perfIdx] == cellIdx) {
960 for (
int i = 0; i < RateVector::dimension; ++i) {
961 rates[i] += connectionRates_[perfIdx][i];
967 template<
typename TypeTag>
971 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
972 if (this->cells()[perfIdx] == cellIdx) {
973 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
974 return connectionRates_[perfIdx][activeCompIdx].value();
978 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
986 template<
typename TypeTag>
994 if (!this->param_.check_well_operability_) {
998 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
1002 updateWellOperability(simulator, well_state, deferred_logger);
1003 if (!this->operability_status_.isOperableAndSolvable()) {
1004 this->operability_status_.use_vfpexplicit =
true;
1005 deferred_logger.
debug(
"EXPLICIT_LOOKUP_VFP",
1006 "well not operable, trying with explicit vfp lookup: " + this->name());
1007 updateWellOperability(simulator, well_state, deferred_logger);
1013 template<
typename TypeTag>
1024 const auto& summary_state = simulator.vanguard().summaryState();
1025 const auto& well_name = this->name();
1026 if (!this->wellHasTHPConstraints(summary_state)) {
1027 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1028 deferred_logger.
info(msg);
1031 const auto& schedule = simulator.vanguard().schedule();
1032 const auto report_step_idx = simulator.episodeIndex();
1033 const auto& glo = schedule.glo(report_step_idx);
1034 if (!glo.has_well(well_name)) {
1035 const std::string msg = fmt::format(
1036 "GLIFT WTEST: Well {} : Gas lift not activated: "
1037 "WLIFTOPT is probably missing. Skipping.", well_name);
1038 deferred_logger.
info(msg);
1041 const auto& gl_well = glo.well(well_name);
1044 std::unique_ptr<GasLiftSingleWell> glift =
1045 initializeGliftWellTest_<GasLiftSingleWell>(simulator,
1051 auto [wtest_alq, success] = glift->wellTestALQ();
1053 const auto& unit_system = schedule.getUnits();
1055 well_state.
well(well_name).alq_state.set(wtest_alq);
1057 "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1058 well_name, unit_system.from_si(UnitSystem::measure::gas_surface_rate, wtest_alq));
1061 if (!gl_well.use_glo()) {
1063 "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1065 unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.
well(well_name).alq_state.get()));
1070 "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1074 deferred_logger.
info(msg);
1077 template<
typename TypeTag>
1085 if (this->param_.local_well_solver_control_switching_) {
1086 const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger);
1090 deferred_logger.
debug(
"Operability check using well equations did not converge for well "
1091 + this->name() +
", reverting to classical approach." );
1094 this->operability_status_.resetOperability();
1096 bool thp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1097 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1098 bool bhp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1099 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1103 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
1104 if (check_thp || bhp_controlled) {
1105 updateIPR(simulator, deferred_logger);
1106 checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
1110 checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger);
1114 template<
typename TypeTag>
1123 assert(this->param_.local_well_solver_control_switching_);
1124 this->operability_status_.resetOperability();
1126 const auto& group_state = simulator.problem().wellModel().groupState();
1127 const double dt = simulator.timeStepSize();
1129 bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger);
1133 template<
typename TypeTag>
1143 const auto& well = this->well_ecl_;
1144 const int well_index = this->index_of_well_;
1145 auto& ws = well_state.
well(well_index);
1148 const auto& summaryState = simulator.vanguard().summaryState();
1149 const auto& schedule = simulator.vanguard().schedule();
1153 ws.primaryvar.resize(0);
1155 if (this->wellIsStopped()) {
1156 for (
int p = 0; p<np; ++p) {
1157 ws.surface_rates[p] = 0;
1163 if (this->isInjector() )
1165 const auto& controls = well.injectionControls(summaryState);
1167 InjectorType injectorType = controls.injector_type;
1169 switch (injectorType) {
1170 case InjectorType::WATER:
1175 case InjectorType::OIL:
1180 case InjectorType::GAS:
1186 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1189 const auto current = ws.injection_cmode;
1192 case Well::InjectorCMode::RATE:
1194 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1195 if(this->rsRvInj() > 0) {
1196 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1198 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1201 OPM_DEFLOG_THROW(std::runtime_error,
"Expected OIL or GAS as type for injectors when RS/RV (item 10) is non-zero " + this->name(), deferred_logger );
1207 case Well::InjectorCMode::RESV:
1209 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1210 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
1211 const Scalar coeff = convert_coeff[phasePos];
1212 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1216 case Well::InjectorCMode::THP:
1218 auto rates = ws.surface_rates;
1223 this->getRefDensity(),
1226 ws.thp = this->getTHPConstraint(summaryState);
1231 Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1232 if (total_rate <= 0.0)
1233 ws.surface_rates = ws.well_potentials;
1237 case Well::InjectorCMode::BHP:
1239 ws.bhp = controls.bhp_limit;
1241 for (
int p = 0; p<np; ++p) {
1242 total_rate += ws.surface_rates[p];
1247 if (total_rate <= 0.0)
1248 ws.surface_rates = ws.well_potentials;
1252 case Well::InjectorCMode::GRUP:
1254 assert(well.isAvailableForGroupControl());
1255 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1256 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1257 well_state[well.name()].efficiency_scaling_factor;
1258 std::optional<Scalar> target =
1259 this->getGroupInjectionTargetRate(group,
1268 ws.surface_rates[phasePos] = *target;
1271 case Well::InjectorCMode::CMODE_UNDEFINED:
1273 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
1283 ws.surface_rates[phasePos] = std::max(
Scalar{1.e-7}, ws.surface_rates[phasePos]);
1286 ws.bhp = controls.bhp_limit;
1292 const auto current = ws.production_cmode;
1293 const auto& controls = well.productionControls(summaryState);
1295 case Well::ProducerCMode::ORAT:
1297 Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
1300 if (current_rate > 0.0) {
1301 for (
int p = 0; p<np; ++p) {
1302 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1305 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1306 double control_fraction = fractions[pu.phase_pos[Oil]];
1307 if (control_fraction != 0.0) {
1308 for (
int p = 0; p<np; ++p) {
1309 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1315 case Well::ProducerCMode::WRAT:
1317 Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
1320 if (current_rate > 0.0) {
1321 for (
int p = 0; p<np; ++p) {
1322 ws.surface_rates[p] *= controls.water_rate/current_rate;
1325 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1326 const Scalar control_fraction = fractions[pu.phase_pos[Water]];
1327 if (control_fraction != 0.0) {
1328 for (
int p = 0; p<np; ++p) {
1329 ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
1335 case Well::ProducerCMode::GRAT:
1337 Scalar current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
1340 if (current_rate > 0.0) {
1341 for (
int p = 0; p<np; ++p) {
1342 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1345 const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
1346 const Scalar control_fraction = fractions[pu.phase_pos[Gas]];
1347 if (control_fraction != 0.0) {
1348 for (
int p = 0; p<np; ++p) {
1349 ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
1357 case Well::ProducerCMode::LRAT:
1359 Scalar current_rate = - ws.surface_rates[ pu.phase_pos[Water] ]
1360 - ws.surface_rates[ pu.phase_pos[Oil] ];
1363 if (current_rate > 0.0) {
1364 for (
int p = 0; p<np; ++p) {
1365 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1368 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1369 const Scalar control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
1370 if (control_fraction != 0.0) {
1371 for (
int p = 0; p<np; ++p) {
1372 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1378 case Well::ProducerCMode::CRAT:
1381 fmt::format(
"CRAT control not supported, well {}", this->name()),
1384 case Well::ProducerCMode::RESV:
1386 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1387 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1388 Scalar total_res_rate = 0.0;
1389 for (
int p = 0; p<np; ++p) {
1390 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1392 if (controls.prediction_mode) {
1395 if (total_res_rate > 0.0) {
1396 for (
int p = 0; p<np; ++p) {
1397 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1400 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1401 for (
int p = 0; p<np; ++p) {
1402 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1406 std::vector<Scalar> hrates(this->number_of_phases_,0.);
1407 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1408 hrates[pu.phase_pos[Water]] = controls.water_rate;
1410 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1411 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
1413 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1414 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
1416 std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
1417 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1418 Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1421 if (total_res_rate > 0.0) {
1422 for (
int p = 0; p<np; ++p) {
1423 ws.surface_rates[p] *= target/total_res_rate;
1426 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1427 for (
int p = 0; p<np; ++p) {
1428 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1434 case Well::ProducerCMode::BHP:
1436 ws.bhp = controls.bhp_limit;
1438 for (
int p = 0; p<np; ++p) {
1439 total_rate -= ws.surface_rates[p];
1444 if (total_rate <= 0.0){
1445 for (
int p = 0; p<np; ++p) {
1446 ws.surface_rates[p] = -ws.well_potentials[p];
1451 case Well::ProducerCMode::THP:
1453 const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger);
1455 if (!update_success) {
1459 auto rates = ws.surface_rates;
1460 this->adaptRatesForVFP(rates);
1462 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1464 ws.thp = this->getTHPConstraint(summaryState);
1468 const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1469 if (total_rate <= 0.0) {
1470 for (
int p = 0; p < this->number_of_phases_; ++p) {
1471 ws.surface_rates[p] = -ws.well_potentials[p];
1477 case Well::ProducerCMode::GRUP:
1479 assert(well.isAvailableForGroupControl());
1480 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1481 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1482 well_state[well.name()].efficiency_scaling_factor;
1483 Scalar scale = this->getGroupProductionTargetRate(group,
1493 for (
int p = 0; p<np; ++p) {
1494 ws.surface_rates[p] *= scale;
1496 ws.trivial_group_target =
false;
1500 ws.trivial_group_target =
true;
1504 case Well::ProducerCMode::CMODE_UNDEFINED:
1507 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1513 ws.bhp = controls.bhp_limit;
1518 template<
typename TypeTag>
1527 const bool isGroupControlled = this->wellUnderGroupControl(well_state.
well(this->index_of_well_));
1528 if (!isGroupControlled) {
1530 const auto& summaryState = simulator.vanguard().summaryState();
1531 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1533 return this->wellUnderZeroGroupRateTarget(simulator, well_state, deferred_logger, isGroupControlled);
1537 template <
typename TypeTag>
1542 const std::optional<bool> group_control)
const
1545 const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.
well(this->index_of_well_)));
1546 if (isGroupControlled) {
1547 const auto& summaryState = simulator.vanguard().summaryState();
1548 const auto& group_state = simulator.problem().wellModel().groupState();
1549 const auto& schedule = simulator.vanguard().schedule();
1550 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state, deferred_logger);
1555 template<
typename TypeTag>
1564 return this->wellIsStopped()
1565 || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
1568 template<
typename TypeTag>
1569 std::vector<typename WellInterface<TypeTag>::Scalar>
1575 const int np = this->number_of_phases_;
1576 std::vector<Scalar> scaling_factor(np);
1577 const auto& ws = well_state.
well(this->index_of_well_);
1579 Scalar total_potentials = 0.0;
1580 for (
int p = 0; p<np; ++p) {
1581 total_potentials += ws.well_potentials[p];
1583 if (total_potentials > 0) {
1584 for (
int p = 0; p<np; ++p) {
1585 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1587 return scaling_factor;
1592 const int nperf = this->number_of_local_perforations_;
1593 for (
int perf = 0; perf < nperf; ++perf) {
1594 total_tw += this->well_index_[perf];
1596 total_tw = this->parallelWellInfo().communication().sum(total_tw);
1598 for (
int perf = 0; perf < nperf; ++perf) {
1599 const int cell_idx = this->well_cells_[perf];
1600 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1601 const auto& fs = intQuants.fluidState();
1602 const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
1603 Scalar total_mobility = 0.0;
1604 for (
int p = 0; p < np; ++p) {
1605 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1606 total_mobility += fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value();
1608 for (
int p = 0; p < np; ++p) {
1609 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1610 scaling_factor[p] += well_tw_fraction * fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value() / total_mobility;
1613 return scaling_factor;
1618 template <
typename TypeTag>
1625 assert(this->isProducer());
1629 auto& ws = well_state.
well(this->index_of_well_);
1630 int nonzero_rate_index = -1;
1631 const Scalar floating_point_error_epsilon = 1e-14;
1632 for (
int p = 0; p < this->number_of_phases_; ++p) {
1633 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1634 if (nonzero_rate_index == -1) {
1635 nonzero_rate_index = p;
1644 std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
1645 bool rates_evaluated_at_1bar =
false;
1647 const auto& summary_state = simulator.vanguard().summaryState();
1648 const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
1649 const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
1650 this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
1652 rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
1654 if (std::any_of(well_q_s.begin(), well_q_s.end(), [](
Scalar q) { return q > 0.0; })) {
1656 if (!rates_evaluated_at_1bar) {
1657 this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
1658 rates_evaluated_at_1bar =
true;
1661 for (
auto& q : well_q_s) {
1662 q = std::min(q,
Scalar{0.0});
1667 if (nonzero_rate_index == -1) {
1671 const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
1672 for (
int p = 0; p < this->number_of_phases_; ++p) {
1673 ws.surface_rates[p] = factor * well_q_s[p];
1682 const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1683 const Scalar computed_rate = well_q_s[nonzero_rate_index];
1684 if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {
1686 const Scalar factor = initial_nonzero_rate / computed_rate;
1687 assert(factor < 1.0);
1688 for (
int p = 0; p < this->number_of_phases_; ++p) {
1690 if (p != nonzero_rate_index) {
1691 ws.surface_rates[p] = factor * well_q_s[p];
1700 for (
int p = 0; p < this->number_of_phases_; ++p) {
1701 ws.surface_rates[p] = well_q_s[p];
1705 template <
typename TypeTag>
1706 std::vector<typename WellInterface<TypeTag>::Scalar>
1713 OPM_TIMEFUNCTION_LOCAL();
1716 if (
static_cast<std::size_t
>(perf) >= this->well_cells_.size()) {
1717 OPM_THROW(std::invalid_argument,
"The perforation index exceeds the size of the local containers - possibly wellIndex was called with a global instead of a local perforation index!");
1719 auto wi = std::vector<Scalar>
1720 (this->num_components_, this->well_index_[perf] * trans_mult);
1722 if constexpr (! Indices::gasEnabled) {
1726 const auto& wdfac = this->well_ecl_.getWDFAC();
1728 if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1732 const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
1739 const auto& connection = this->well_ecl_.getConnections()[ws.
perf_data.ecl_index[perf]];
1740 const Scalar Kh = connection.Kh();
1741 const Scalar scaling = 3.141592653589 * Kh * connection.wpimult();
1742 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1745 const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1746 const Scalar drawdown = cell_pressure - connection_pressure;
1747 const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1748 const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1750 const Scalar b = 2*scaling/wi[gas_comp_idx];
1751 const Scalar c = -2*scaling*mob_g*drawdown;
1753 Scalar consistent_Q = -1.0e20;
1755 const Scalar r2n = b*b + 4*a*c;
1757 const Scalar rn = std::sqrt(r2n);
1758 const Scalar xn1 = (b-rn)*0.5/a;
1762 const Scalar xn2 = (b+rn)*0.5/a;
1763 if (xn2 <= 0 && xn2 > consistent_Q) {
1769 const Scalar r2p = b*b - 4*a*c;
1771 const Scalar rp = std::sqrt(r2p);
1772 const Scalar xp1 = (rp-b)*0.5/a;
1773 if (xp1 > 0 && xp1 < consistent_Q) {
1776 const Scalar xp2 = -(rp+b)*0.5/a;
1777 if (xp2 > 0 && xp2 < consistent_Q) {
1781 wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1786 template <
typename TypeTag>
1792 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1796 auto& d_factor = ws.
perf_data.connection_d_factor;
1798 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1799 const int cell_idx = this->well_cells_[perf];
1800 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1802 d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
1806 template <
typename TypeTag>
1813 auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
1814 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
1818 auto gas_visc = [connection_pressure = ws.
perf_data.pressure[perf],
1820 regIdx = this->pvtRegionIdx(), &intQuants]()
1822 const auto rv = getValue(intQuants.fluidState().Rv());
1824 const auto& gasPvt = FluidSystem::gasPvt();
1829 const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
1830 (regIdx, temperature, connection_pressure);
1832 if (! (rv < rv_sat)) {
1833 return gasPvt.saturatedViscosity(regIdx, temperature,
1834 connection_pressure);
1837 return gasPvt.viscosity(regIdx, temperature, connection_pressure,
1838 rv, getValue(intQuants.fluidState().Rvw()));
1841 const auto& connection = this->well_ecl_.getConnections()
1844 return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
1848 template <
typename TypeTag>
1854 auto connCF = [&connIx = std::as_const(ws.
perf_data.ecl_index),
1855 &conns = this->well_ecl_.getConnections()]
1858 return conns[connIx[perf]].CF();
1861 auto& tmult = ws.
perf_data.connection_compaction_tmult;
1862 auto& ctf = ws.
perf_data.connection_transmissibility_factor;
1864 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1865 const int cell_idx = this->well_cells_[perf];
1867 const auto& intQuants = simulator.model()
1868 .intensiveQuantities(cell_idx, 0);
1870 tmult[perf] = simulator.problem()
1871 .template wellTransMultiplier<double>(intQuants, cell_idx);
1873 ctf[perf] = connCF(perf) * tmult[perf];
1878 template<
typename TypeTag>
1882 if constexpr (Indices::oilEnabled) {
1883 return fs.pressure(FluidSystem::oilPhaseIdx);
1884 }
else if constexpr (Indices::gasEnabled) {
1885 return fs.pressure(FluidSystem::gasPhaseIdx);
1887 return fs.pressure(FluidSystem::waterPhaseIdx);
1891 template <
typename TypeTag>
1892 template<
class Value,
class Callback>
1896 const int local_perf_index,
1897 std::vector<Value>& mob,
1898 Callback& extendEval,
1901 auto relpermArray = []()
1903 if constexpr (std::is_same_v<Value, Scalar>) {
1904 return std::array<Scalar,3>{};
1906 return std::array<Eval,3>{};
1909 if (
static_cast<std::size_t
>(local_perf_index) >= this->well_cells_.size()) {
1910 OPM_THROW(std::invalid_argument,
"The perforation index exceeds the size of the local containers - possibly getMobility was called with a global instead of a local perforation index!");
1912 const int cell_idx = this->well_cells_[local_perf_index];
1913 assert (
int(mob.size()) == this->num_components_);
1914 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1915 const auto& materialLawManager = simulator.problem().materialLawManager();
1919 const int satid = this->saturation_table_number_[local_perf_index] - 1;
1920 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1921 if (satid == satid_elem) {
1922 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1923 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1927 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1928 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1930 if constexpr (has_solvent) {
1931 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1934 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1935 auto relativePerms = relpermArray();
1936 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1939 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1942 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1943 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1947 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1948 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1952 if constexpr (has_solvent) {
1953 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1957 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1958 const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
1959 const auto& connections = this->well_ecl_.getConnections();
1960 const auto& connection = connections[perf_ecl_index];
1961 if (connection.filterCakeActive()) {
1962 std::transform(mob.begin(), mob.end(), mob.begin(),
1963 [mult = this->inj_fc_multiplier_[local_perf_index] ](
const auto val)
1964 { return val * mult; });
1970 template<
typename TypeTag>
1978 const auto& summary_state = simulator.vanguard().summaryState();
1980 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1981 simulator, summary_state, this->getALQ(well_state), deferred_logger,
false);
1982 if (bhp_at_thp_limit) {
1983 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
1984 if (thp_update_iterations) {
1985 computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
1986 rates, deferred_logger);
1988 computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
1989 rates, deferred_logger);
1991 auto& ws = well_state.
well(this->name());
1992 ws.surface_rates = rates;
1993 ws.bhp = *bhp_at_thp_limit;
1994 ws.thp = this->getTHPConstraint(summary_state);
2001 template <
typename TypeTag>
2006 const std::vector<Scalar>& mobility,
2010 const int np = this->number_of_phases_;
2011 for (
int p = 0; p < np; ++p) {
2014 const auto connMob =
2015 mobility[this->flowPhaseToModelCompIdx(p)]
2016 * fs.invB(this->flowPhaseToModelPhaseIdx(p)).value();
2018 connPI[p] = connPICalc(connMob);
2021 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2022 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2024 const auto io = pu.phase_pos[Oil];
2025 const auto ig = pu.phase_pos[Gas];
2027 const auto vapoil = connPI[ig] * fs.Rv().value();
2028 const auto disgas = connPI[io] * fs.Rs().value();
2030 connPI[io] += vapoil;
2031 connPI[ig] += disgas;
2036 template <
typename TypeTag>
2040 const Phase preferred_phase,
2042 const std::vector<Scalar>& mobility,
2050 if (preferred_phase == Phase::GAS) {
2051 phase_pos = pu.phase_pos[Gas];
2053 else if (preferred_phase == Phase::OIL) {
2054 phase_pos = pu.phase_pos[Oil];
2056 else if (preferred_phase == Phase::WATER) {
2057 phase_pos = pu.phase_pos[Water];
2061 fmt::format(
"Unsupported Injector Type ({}) "
2062 "for well {} during connection I.I. calculation",
2063 static_cast<int>(preferred_phase), this->name()),
2067 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2068 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToModelPhaseIdx(phase_pos)).value());
2071 template<
typename TypeTag>
2072 template<
class GasLiftSingleWell>
2073 std::unique_ptr<GasLiftSingleWell>
2083 auto& comm = simulator.vanguard().grid().comm();
2084 ecl_well_map.try_emplace(this->name(), &(this->wellEcl()), this->indexOfWell());
2087 simulator.vanguard().schedule(),
2088 simulator.vanguard().summaryState(),
2089 simulator.episodeIndex(),
2090 simulator.model().newtonMethod().numIterations(),
2100 std::set<int> sync_groups;
2101 const auto& summary_state = simulator.vanguard().summaryState();
2102 return std::make_unique<GasLiftSingleWell>(*
this,
#define OPM_DEFLOG_THROW(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:45
@ Liquid
Definition: BlackoilPhases.hpp:42
@ Aqua
Definition: BlackoilPhases.hpp:42
@ Vapour
Definition: BlackoilPhases.hpp:42
Definition: DeferredLogger.hpp:57
void info(const std::string &tag, const std::string &message)
void warning(const std::string &tag, const std::string &message)
void debug(const std::string &tag, const std::string &message)
Definition: GasLiftGroupInfo.hpp:46
Definition: GroupState.hpp:43
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:195
Definition: SingleWellState.hpp:42
PerfData< Scalar > perf_data
Definition: SingleWellState.hpp:115
Scalar temperature
Definition: SingleWellState.hpp:94
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
std::optional< Scalar > estimateStableBhp(const WellState< Scalar > &well_state, const Well &well, const std::vector< Scalar > &rates, const Scalar rho, const SummaryState &summaryState) const
Scalar calculateMinimumBhpFromThp(const WellState< Scalar > &well_state, const Well &well, const SummaryState &summaryState, const Scalar rho) const
bool isStableSolution(const WellState< Scalar > &well_state, const Well &well, const std::vector< Scalar > &rates, const SummaryState &summaryState) const
EvalWell calculateBhpFromThp(const WellState< Scalar > &well_state, const std::vector< EvalWell > &rates, const Well &well, const SummaryState &summaryState, const Scalar rho, DeferredLogger &deferred_logger) const
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
int number_of_local_perforations_
Definition: WellInterfaceGeneric.hpp:331
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:297
const std::vector< FluidSystem::Scalar > & wellIndex() const
Definition: WellInterfaceGeneric.hpp:146
FluidSystem::Scalar wsolvent_
Definition: WellInterfaceGeneric.hpp:373
Definition: WellInterfaceIndices.hpp:34
bool updateWellControl(const Simulator &simulator, const IndividualOrGroup iog, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:191
Scalar woxygen() const
Definition: WellInterface_impl.hpp:167
IndividualOrGroup
Definition: WellInterface.hpp:242
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:82
WellInterface(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters ¶m, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Constructor.
Definition: WellInterface_impl.hpp:58
void assembleWellEq(const Simulator &simulator, const double dt, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:841
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:100
BlackOilFluidState< Eval, FluidSystem, has_temperature, has_energy, Indices::compositionSwitchIdx >=0, has_watVapor, has_brine, has_saltPrecip, has_disgas_in_water, Indices::numPhases > FluidState
Definition: WellInterface.hpp:132
Scalar wfoam() const
Definition: WellInterface_impl.hpp:129
std::vector< RateVector > connectionRates_
Definition: WellInterface.hpp:367
bool wellUnderZeroRateTarget(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1521
void computeConnLevelProdInd(const FluidState &fs, const std::function< Scalar(const Scalar)> &connPICalc, const std::vector< Scalar > &mobility, Scalar *connPI) const
Definition: WellInterface_impl.hpp:2004
bool updateWellStateWithTHPTargetProd(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1973
void addCellRates(RateVector &rates, int cellIdx) const
Definition: WellInterface_impl.hpp:953
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const
Definition: WellInterface_impl.hpp:969
void getMobility(const Simulator &simulator, const int local_perf_index, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1895
bool iterateWellEquations(const Simulator &simulator, const double dt, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:503
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:86
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:83
void checkWellOperability(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:989
std::optional< Scalar > estimateOperableBhp(const Simulator &ebos_simulator, const double dt, WellState< Scalar > &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:666
void updateWellOperability(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1080
void wellTesting(const Simulator &simulator, const double simulation_time, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, WellTestState &welltest_state, const PhaseUsage &phase_usage, GLiftEclWells &ecl_well_map, std::map< std::string, double > &open_times, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:381
bool solveWellWithOperabilityCheck(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:566
virtual void updateWellStateWithTarget(const Simulator &simulator, const GroupState< Scalar > &group_state, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1136
Eval getPerfCellPressure(const FluidState &fs) const
Definition: WellInterface_impl.hpp:1880
void updateConnectionTransmissibilityFactor(const Simulator &simulator, SingleWellState< Scalar > &ws) const
Definition: WellInterface_impl.hpp:1851
bool stoppedOrZeroRateTarget(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1558
void initializeProducerWellState(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1621
typename Base::ModelParameters ModelParameters
Definition: WellInterface.hpp:106
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:84
static constexpr bool has_solvent
Definition: WellInterface.hpp:108
void prepareWellBeforeAssembling(const Simulator &simulator, const double dt, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:877
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: WellInterface.hpp:89
void computeConnLevelInjInd(const FluidState &fs, const Phase preferred_phase, const std::function< Scalar(const Scalar)> &connIICalc, const std::vector< Scalar > &mobility, Scalar *connII, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:2039
void solveWellEquation(const Simulator &simulator, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:779
virtual void init(const PhaseUsage *phase_usage_arg, const std::vector< Scalar > &depth_arg, const Scalar gravity_arg, const std::vector< Scalar > &B_avg, const bool changed_to_open_this_step)
Definition: WellInterface_impl.hpp:95
void updateConnectionDFactor(const Simulator &simulator, SingleWellState< Scalar > &ws) const
Definition: WellInterface_impl.hpp:1789
std::vector< Scalar > initialWellRateFractions(const Simulator &ebosSimulator, const WellState< Scalar > &well_state) const
Definition: WellInterface_impl.hpp:1571
Scalar wsalt() const
Definition: WellInterface_impl.hpp:143
bool updateWellOperabilityFromWellEq(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1117
void gliftBeginTimeStepWellTestUpdateALQ(const Simulator &simulator, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, const PhaseUsage &phase_usage, GLiftEclWells &ecl_well_map, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1016
typename GasLiftGroupInfo< Scalar >::GLiftEclWells GLiftEclWells
Definition: WellInterface.hpp:91
typename Base::Eval Eval
Definition: WellInterface.hpp:95
Scalar wpolymer() const
Definition: WellInterface_impl.hpp:113
Scalar computeConnectionDFactor(const int perf, const IntensiveQuantities &intQuants, const SingleWellState< Scalar > &ws) const
Definition: WellInterface_impl.hpp:1809
bool solveWellWithBhp(const Simulator &simulator, const double dt, const Scalar bhp, WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:699
bool solveWellWithZeroRate(const Simulator &simulator, const double dt, WellState< Scalar > &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:736
bool solveWellForTesting(const Simulator &simulator, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:757
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:85
Scalar wurea() const
Definition: WellInterface_impl.hpp:179
bool wellUnderZeroGroupRateTarget(const Simulator &simulator, const WellState< Scalar > &well_state, DeferredLogger &deferred_logger, std::optional< bool > group_control=std::nullopt) const
Definition: WellInterface_impl.hpp:1539
Scalar wmicrobes() const
Definition: WellInterface_impl.hpp:155
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:857
bool updateWellControlAndStatusLocalIteration(const Simulator &simulator, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, const Scalar WQTotal, DeferredLogger &deferred_logger, const bool fixed_control=false, const bool fixed_status=false)
Definition: WellInterface_impl.hpp:274
std::unique_ptr< GasLiftSingleWell > initializeGliftWellTest_(const Simulator &simulator, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, const PhaseUsage &phase_usage, GLiftEclWells &ecl_well_map, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:2075
static constexpr bool has_zFraction
Definition: WellInterface.hpp:109
Definition: WellState.hpp:65
const SingleWellState< Scalar > & well(std::size_t well_index) const
Definition: WellState.hpp:285
int numPhases() const
The number of phases present.
Definition: WellState.hpp:239
@ NONE
Definition: DeferredLogger.hpp:46
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:39
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Static data associated with a well perforation.
Definition: PerforationData.hpp:30
Definition: BlackoilPhases.hpp:46