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_conservation_quantities,
66 const int index_of_well,
74 num_conservation_quantities,
82 if (well.isInjector()) {
83 auto injectorType = this->
well_ecl_.injectorType();
84 if (injectorType == InjectorType::GAS) {
92 template<
typename TypeTag>
95 init(
const std::vector<Scalar>& ,
97 const std::vector<Scalar>& B_avg,
98 const bool changed_to_open_this_step)
100 this->gravity_ = gravity_arg;
102 this->changed_to_open_this_step_ = changed_to_open_this_step;
108 template<
typename TypeTag>
113 if constexpr (has_polymer) {
114 return this->wpolymer_();
124 template<
typename TypeTag>
129 if constexpr (has_foam) {
130 return this->wfoam_();
138 template<
typename TypeTag>
143 if constexpr (has_brine) {
144 return this->wsalt_();
150 template<
typename TypeTag>
155 if constexpr (has_micp) {
156 return this->wmicrobes_();
162 template<
typename TypeTag>
167 if constexpr (has_micp) {
168 return this->woxygen_();
174 template<
typename TypeTag>
179 if constexpr (has_micp) {
180 return this->wurea_();
186 template<
typename TypeTag>
196 if (stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
200 const auto& summaryState = simulator.vanguard().summaryState();
201 const auto& schedule = simulator.vanguard().schedule();
202 const auto& well = this->well_ecl_;
203 auto& ws = well_state.
well(this->index_of_well_);
205 if (well.isInjector()) {
206 from = WellInjectorCMode2String(ws.injection_cmode);
208 from = WellProducerCMode2String(ws.production_cmode);
211 const int episodeIdx = simulator.episodeIndex();
212 const int iterationIdx = simulator.model().newtonMethod().numIterations();
213 const int nupcol = schedule[episodeIdx].nupcol();
214 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
217 const bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
219 const auto msg = fmt::format(
" The control mode for well {} is oscillating. \n"
220 "We don't allow for more than {} switches after NUPCOL iterations. (NUPCOL = {}) \n"
221 "The control is kept at {}.",
222 this->name(), this->param_.max_number_of_well_switches_, nupcol, from);
223 deferred_logger.
info(msg);
225 this->well_control_log_.push_back(from);
229 bool changed =
false;
230 if (iog == IndividualOrGroup::Individual) {
231 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
232 }
else if (iog == IndividualOrGroup::Group) {
233 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState,
true, deferred_logger);
235 assert(iog == IndividualOrGroup::Both);
236 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
242 if (well.isInjector()) {
243 to = WellInjectorCMode2String(ws.injection_cmode);
245 to = WellProducerCMode2String(ws.production_cmode);
247 std::ostringstream ss;
248 ss <<
" Switching control mode for well " << this->name()
252 ss <<
" on rank " << cc.rank();
254 deferred_logger.
debug(ss.str());
259 if (iterationIdx >= nupcol || this->well_control_log_.empty()) {
260 this->well_control_log_.push_back(from);
262 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
263 updatePrimaryVariables(simulator, well_state, deferred_logger);
269 template<
typename TypeTag>
275 const Well::InjectionControls& inj_controls,
276 const Well::ProductionControls& prod_controls,
279 const bool fixed_control,
280 const bool fixed_status)
283 const auto& summary_state = simulator.vanguard().summaryState();
284 const auto& schedule = simulator.vanguard().schedule();
285 auto& ws = well_state.
well(this->index_of_well_);
287 if (this->isInjector()) {
288 from = WellInjectorCMode2String(ws.injection_cmode);
290 from = WellProducerCMode2String(ws.production_cmode);
292 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
294 if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(well_state.
well(this->index_of_well_).status == WellStatus::OPEN)) {
298 const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
299 if (!this->wellIsStopped()){
300 if (wqTotal*sgn <= 0.0 && !fixed_status){
304 bool changed =
false;
305 if (!fixed_control) {
310 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
311 prod_controls.hasControl(Well::ProducerCMode::GRUP);
312 bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
313 if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
314 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
316 if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
317 changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,
false, deferred_logger);
321 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
322 ws.production_cmode == Well::ProducerCMode::THP;
324 ws.thp = this->getTHPConstraint(summary_state);
327 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
329 updatePrimaryVariables(simulator, well_state, deferred_logger);
334 }
else if (!fixed_status){
336 const Scalar bhp = well_state.
well(this->index_of_well_).bhp;
337 Scalar prod_limit = prod_controls.bhp_limit;
338 Scalar inj_limit = inj_controls.bhp_limit;
339 const bool has_thp = this->wellHasTHPConstraints(summary_state);
341 std::vector<Scalar> rates(this->num_conservation_quantities_);
342 if (this->isInjector()){
344 calculateBhpFromThp(well_state, rates,
347 this->getRefDensity(),
349 inj_limit = std::min(bhp_thp,
static_cast<Scalar>(inj_controls.bhp_limit));
354 calculateMinimumBhpFromThp(well_state,
357 this->getRefDensity());
358 prod_limit = std::max(bhp_min,
static_cast<Scalar>(prod_controls.bhp_limit));
361 const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
364 well_state.
well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
366 well_state.
well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
377 template<
typename TypeTag>
381 const double simulation_time,
384 WellTestState& well_test_state,
386 std::map<std::string, double>& open_times,
390 deferred_logger.
info(
" well " + this->name() +
" is being tested");
393 auto& ws = well_state_copy.
well(this->indexOfWell());
395 const auto& summary_state = simulator.vanguard().summaryState();
396 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
397 if (this->isProducer()) {
398 ws.production_cmode = has_thp_limit ? Well::ProducerCMode::THP : Well::ProducerCMode::BHP;
400 ws.injection_cmode = has_thp_limit ? Well::InjectorCMode::THP : Well::InjectorCMode::BHP;
405 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
406 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
407 updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
409 if (this->isProducer()) {
410 const auto& schedule = simulator.vanguard().schedule();
411 const auto report_step = simulator.episodeIndex();
412 const auto& glo = schedule.glo(report_step);
414 gliftBeginTimeStepWellTestUpdateALQ(simulator,
422 WellTestState welltest_state_temp;
424 bool testWell =
true;
429 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
430 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
432 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
433 deferred_logger.
debug(msg);
438 updateWellOperability(simulator, well_state_copy, deferred_logger);
439 if ( !this->isOperableAndSolvable() ) {
440 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
441 deferred_logger.
debug(msg);
444 std::vector<Scalar> potentials;
446 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
447 }
catch (
const std::exception& e) {
448 const std::string msg = fmt::format(
"well {}: computeWellPotentials() "
449 "failed during testing for re-opening: ",
450 this->name(), e.what());
451 deferred_logger.
info(msg);
454 const int np = well_state_copy.
numPhases();
455 for (
int p = 0; p < np; ++p) {
456 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
458 const bool under_zero_target = this->wellUnderZeroGroupRateTarget(simulator, well_state_copy, deferred_logger);
459 this->updateWellTestState(well_state_copy.
well(this->indexOfWell()),
465 this->closeCompletions(welltest_state_temp);
471 if ( welltest_state_temp.num_closed_wells() > 0 ||
472 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
478 if (!welltest_state_temp.well_is_closed(this->name())) {
479 well_test_state.open_well(this->name());
481 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
482 deferred_logger.
info(msg);
485 for (
const auto& completion : this->well_ecl_.getCompletions()) {
486 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
487 well_test_state.open_completion(this->name(), completion.first);
489 well_state = well_state_copy;
490 open_times.try_emplace(this->name(), well_test_state.lastTestTime(this->name()));
497 template<
typename TypeTag>
507 const auto& summary_state = simulator.vanguard().summaryState();
508 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
509 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
510 const auto& ws = well_state.
well(this->indexOfWell());
511 const auto pmode_orig = ws.production_cmode;
512 const auto imode_orig = ws.injection_cmode;
513 bool converged =
false;
516 if (!this->param_.local_well_solver_control_switching_){
517 converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
519 if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.
well(this->index_of_well_).status == WellStatus::OPEN)) {
520 converged = solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
522 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
526 }
catch (NumericalProblem& e ) {
527 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
528 deferred_logger.
warning(
"INNER_ITERATION_FAILED", msg);
533 if (ws.production_cmode != pmode_orig || ws.injection_cmode != imode_orig) {
535 if (this->isInjector()) {
536 from = WellInjectorCMode2String(imode_orig);
537 to = WellInjectorCMode2String(ws.injection_cmode);
539 from = WellProducerCMode2String(pmode_orig);
540 to = WellProducerCMode2String(ws.production_cmode);
542 const auto msg = fmt::format(
" Well {} switched from {} to {} during local solve", this->name(), from, to);
543 deferred_logger.
debug(msg);
544 const int episodeIdx = simulator.episodeIndex();
545 const int iterationIdx = simulator.model().newtonMethod().numIterations();
546 const auto& schedule = simulator.vanguard().schedule();
547 const int nupcol = schedule[episodeIdx].nupcol();
551 if (iterationIdx >= nupcol || this->well_control_log_.empty()) {
552 this->well_control_log_.push_back(from);
560 template<
typename TypeTag>
565 const Well::InjectionControls& inj_controls,
566 const Well::ProductionControls& prod_controls,
572 const auto& summary_state = simulator.vanguard().summaryState();
573 bool is_operable =
true;
574 bool converged =
true;
575 auto& ws = well_state.
well(this->index_of_well_);
577 if (this->wellIsStopped()) {
579 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
580 if (!bhp_target.has_value()) {
582 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
583 deferred_logger.
debug(msg);
586 solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
590 ws.thp = this->getTHPConstraint(summary_state);
591 const Scalar bhp = std::max(bhp_target.value(),
592 static_cast<Scalar>(prod_controls.bhp_limit));
593 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
598 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
601 const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
603 if (converged && !stoppedOrZeroRateTarget(simulator, well_state, deferred_logger) && isThp) {
604 auto rates = well_state.
well(this->index_of_well_).surface_rates;
605 this->adaptRatesForVFP(rates);
606 this->updateIPRImplicit(simulator, well_state, deferred_logger);
610 this->operability_status_.use_vfpexplicit =
true;
613 const Scalar reltol = 1e-3;
614 const Scalar cur_bhp = ws.bhp;
615 if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
616 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
617 deferred_logger.
debug(msg);
618 solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger);
620 ws.thp = this->getTHPConstraint(summary_state);
621 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
628 this->operability_status_.use_vfpexplicit =
true;
630 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
631 if (!bhp_target.has_value()) {
635 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
639 const Scalar bhp = std::max(bhp_target.value(),
640 static_cast<Scalar>(prod_controls.bhp_limit));
641 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
642 ws.thp = this->getTHPConstraint(summary_state);
643 const auto msg = fmt::format(
"Well {} did not converge, re-solving with explicit fractions for VFP caculations.", this->name());
644 deferred_logger.
debug(msg);
645 converged = this->iterateWellEqWithSwitching(simulator, dt,
654 is_operable = is_operable && !this->wellIsStopped();
655 this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
656 this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
660 template<
typename TypeTag>
661 std::optional<typename WellInterface<TypeTag>::Scalar>
666 const SummaryState& summary_state,
669 if (!this->wellHasTHPConstraints(summary_state)) {
671 const bool converged = solveWellWithBhp(simulator, dt, bhp_limit, well_state, deferred_logger);
672 if (!converged || this->wellIsStopped()) {
683 const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
684 if (!converged || this->wellIsStopped()) {
687 this->updateIPRImplicit(simulator, well_state, deferred_logger);
688 auto rates = well_state.
well(this->index_of_well_).surface_rates;
689 this->adaptRatesForVFP(rates);
693 template<
typename TypeTag>
705 auto inj_controls = Well::InjectionControls(0);
706 auto prod_controls = Well::ProductionControls(0);
707 auto& ws = well_state.
well(this->index_of_well_);
708 auto cmode_inj = ws.injection_cmode;
709 auto cmode_prod = ws.production_cmode;
710 if (this->isInjector()) {
711 inj_controls.addControl(Well::InjectorCMode::BHP);
712 inj_controls.bhp_limit = bhp;
713 inj_controls.cmode = Well::InjectorCMode::BHP;
714 ws.injection_cmode = Well::InjectorCMode::BHP;
716 prod_controls.addControl(Well::ProducerCMode::BHP);
717 prod_controls.bhp_limit = bhp;
718 prod_controls.cmode = Well::ProducerCMode::BHP;
719 ws.production_cmode = Well::ProducerCMode::BHP;
724 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true);
725 ws.injection_cmode = cmode_inj;
726 ws.production_cmode = cmode_prod;
730 template<
typename TypeTag>
740 const auto well_status_orig = this->wellStatus_;
744 auto inj_controls = Well::InjectionControls(0);
745 auto prod_controls = Well::ProductionControls(0);
746 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true,
true);
747 this->wellStatus_ = well_status_orig;
751 template<
typename TypeTag>
760 const double dt = simulator.timeStepSize();
761 const bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
763 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" converged");
766 const int max_iter = this->param_.max_welleq_iter_;
767 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
773 template<
typename TypeTag>
782 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
787 const double dt = simulator.timeStepSize();
788 bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
798 auto& ws = well_state.
well(this->indexOfWell());
799 bool thp_control =
false;
800 if (this->well_ecl_.isInjector()) {
801 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
803 ws.injection_cmode = Well::InjectorCMode::BHP;
804 if (this->well_control_log_.empty()) {
805 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
809 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
811 ws.production_cmode = Well::ProducerCMode::BHP;
812 if (this->well_control_log_.empty()) {
813 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
818 const std::string msg = std::string(
"The newly opened well ") + this->name()
819 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
820 deferred_logger.
debug(msg);
821 converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
826 const int max_iter = this->param_.max_welleq_iter_;
827 deferred_logger.
debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
829 well_state = well_state0;
835 template <
typename TypeTag>
845 prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
846 assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
851 template <
typename TypeTag>
861 const auto& summary_state = simulator.vanguard().summaryState();
862 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
863 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
866 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
871 template<
typename TypeTag>
881 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
883 if (this->param_.check_well_operability_iter_)
884 checkWellOperability(simulator, well_state, deferred_logger);
887 const int iteration_idx = simulator.model().newtonMethod().numIterations();
888 if (iteration_idx < this->param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
889 const auto& ws = well_state.
well(this->indexOfWell());
890 const bool nonzero_rate_original =
891 std::any_of(ws.surface_rates.begin(),
892 ws.surface_rates.begin() + well_state.
numPhases(),
893 [](
Scalar rate) { return rate != Scalar(0.0); });
895 this->operability_status_.solvable =
true;
896 bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
899 const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
900 if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
904 this->operability_status_.resetOperability();
906 deferred_logger.
debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
910 if (this->param_.shut_unsolvable_wells_) {
911 this->operability_status_.solvable =
false;
915 if (this->operability_status_.has_negative_potentials) {
916 auto well_state_copy = well_state;
917 std::vector<Scalar> potentials;
919 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
920 }
catch (
const std::exception& e) {
921 const std::string msg = fmt::format(
"well {}: computeWellPotentials() failed "
922 "during attempt to recompute potentials for well: ",
923 this->name(), e.what());
924 deferred_logger.
info(msg);
925 this->operability_status_.has_negative_potentials =
true;
927 auto& ws = well_state.
well(this->indexOfWell());
929 for (
int p = 0; p < np; ++p) {
930 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
933 this->changed_to_open_this_step_ =
false;
934 const bool well_operable = this->operability_status_.isOperableAndSolvable();
936 if (!well_operable && old_well_operable) {
937 deferred_logger.
debug(
" well " + this->name() +
" gets STOPPED during iteration ");
939 changed_to_stopped_this_step_ =
true;
940 }
else if (well_operable && !old_well_operable) {
941 deferred_logger.
debug(
" well " + this->name() +
" gets REVIVED during iteration ");
943 changed_to_stopped_this_step_ =
false;
944 this->changed_to_open_this_step_ =
true;
948 template<
typename TypeTag>
952 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
955 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
956 if (this->cells()[perfIdx] == cellIdx) {
957 for (
int i = 0; i < RateVector::dimension; ++i) {
958 rates[i] += connectionRates_[perfIdx][i];
964 template<
typename TypeTag>
968 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
969 if (this->cells()[perfIdx] == cellIdx) {
970 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
971 return connectionRates_[perfIdx][activeCompIdx].value();
975 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
983 template<
typename TypeTag>
991 if (!this->param_.check_well_operability_) {
995 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
999 updateWellOperability(simulator, well_state, deferred_logger);
1000 if (!this->operability_status_.isOperableAndSolvable()) {
1001 this->operability_status_.use_vfpexplicit =
true;
1002 deferred_logger.
debug(
"EXPLICIT_LOOKUP_VFP",
1003 "well not operable, trying with explicit vfp lookup: " + this->name());
1004 updateWellOperability(simulator, well_state, deferred_logger);
1010 template<
typename TypeTag>
1020 const auto& summary_state = simulator.vanguard().summaryState();
1021 const auto& well_name = this->name();
1022 if (!this->wellHasTHPConstraints(summary_state)) {
1023 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1024 deferred_logger.
info(msg);
1027 const auto& schedule = simulator.vanguard().schedule();
1028 const auto report_step_idx = simulator.episodeIndex();
1029 const auto& glo = schedule.glo(report_step_idx);
1030 if (!glo.has_well(well_name)) {
1031 const std::string msg = fmt::format(
1032 "GLIFT WTEST: Well {} : Gas lift not activated: "
1033 "WLIFTOPT is probably missing. Skipping.", well_name);
1034 deferred_logger.
info(msg);
1037 const auto& gl_well = glo.well(well_name);
1040 std::unique_ptr<GasLiftSingleWell> glift =
1041 initializeGliftWellTest_<GasLiftSingleWell>(simulator,
1046 auto [wtest_alq, success] = glift->wellTestALQ();
1048 const auto& unit_system = schedule.getUnits();
1050 well_state.
well(well_name).alq_state.set(wtest_alq);
1052 "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1053 well_name, unit_system.from_si(UnitSystem::measure::gas_surface_rate, wtest_alq));
1056 if (!gl_well.use_glo()) {
1058 "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1060 unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.
well(well_name).alq_state.get()));
1065 "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1069 deferred_logger.
info(msg);
1072 template<
typename TypeTag>
1080 if (this->param_.local_well_solver_control_switching_) {
1081 const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger);
1085 deferred_logger.
debug(
"Operability check using well equations did not converge for well "
1086 + this->name() +
", reverting to classical approach." );
1089 this->operability_status_.resetOperability();
1091 bool thp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1092 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1093 bool bhp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1094 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1098 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
1099 if (check_thp || bhp_controlled) {
1100 updateIPR(simulator, deferred_logger);
1101 checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
1105 checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger);
1109 template<
typename TypeTag>
1118 assert(this->param_.local_well_solver_control_switching_);
1119 this->operability_status_.resetOperability();
1121 const auto& group_state = simulator.problem().wellModel().groupState();
1122 const double dt = simulator.timeStepSize();
1124 bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger);
1128 template<
typename TypeTag>
1138 const auto& well = this->well_ecl_;
1139 const int well_index = this->index_of_well_;
1140 auto& ws = well_state.
well(well_index);
1142 const auto& summaryState = simulator.vanguard().summaryState();
1143 const auto& schedule = simulator.vanguard().schedule();
1147 ws.primaryvar.resize(0);
1149 if (this->wellIsStopped()) {
1150 for (
int p = 0; p<np; ++p) {
1151 ws.surface_rates[p] = 0;
1157 if (this->isInjector() )
1159 const auto& controls = well.injectionControls(summaryState);
1161 InjectorType injectorType = controls.injector_type;
1163 switch (injectorType) {
1164 case InjectorType::WATER:
1166 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1169 case InjectorType::OIL:
1171 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1174 case InjectorType::GAS:
1176 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1180 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1183 const auto current = ws.injection_cmode;
1186 case Well::InjectorCMode::RATE:
1188 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1189 if(this->rsRvInj() > 0) {
1190 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1191 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1192 ws.surface_rates[gas_pos] = controls.surface_rate * this->rsRvInj();
1193 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1194 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1195 ws.surface_rates[oil_pos] = controls.surface_rate * this->rsRvInj();
1197 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 );
1203 case Well::InjectorCMode::RESV:
1205 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1206 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
1207 const Scalar coeff = convert_coeff[phasePos];
1208 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1212 case Well::InjectorCMode::THP:
1214 auto rates = ws.surface_rates;
1219 this->getRefDensity(),
1222 ws.thp = this->getTHPConstraint(summaryState);
1227 Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1228 if (total_rate <= 0.0)
1229 ws.surface_rates = ws.well_potentials;
1233 case Well::InjectorCMode::BHP:
1235 ws.bhp = controls.bhp_limit;
1237 for (
int p = 0; p<np; ++p) {
1238 total_rate += ws.surface_rates[p];
1243 if (total_rate <= 0.0)
1244 ws.surface_rates = ws.well_potentials;
1248 case Well::InjectorCMode::GRUP:
1250 assert(well.isAvailableForGroupControl());
1251 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1252 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1253 well_state[well.name()].efficiency_scaling_factor;
1254 std::optional<Scalar> target =
1255 this->getGroupInjectionTargetRate(group,
1264 ws.surface_rates[phasePos] = *target;
1267 case Well::InjectorCMode::CMODE_UNDEFINED:
1269 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
1279 ws.surface_rates[phasePos] = std::max(
Scalar{1.e-7}, ws.surface_rates[phasePos]);
1282 ws.bhp = controls.bhp_limit;
1288 const auto current = ws.production_cmode;
1289 const auto& controls = well.productionControls(summaryState);
1291 case Well::ProducerCMode::ORAT:
1293 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1294 Scalar current_rate = -ws.surface_rates[oil_pos];
1297 if (current_rate > 0.0) {
1298 for (
int p = 0; p<np; ++p) {
1299 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1302 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1303 double control_fraction = fractions[oil_pos];
1304 if (control_fraction != 0.0) {
1305 for (
int p = 0; p<np; ++p) {
1306 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1312 case Well::ProducerCMode::WRAT:
1314 const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1315 Scalar current_rate = -ws.surface_rates[water_pos];
1318 if (current_rate > 0.0) {
1319 for (
int p = 0; p<np; ++p) {
1320 ws.surface_rates[p] *= controls.water_rate/current_rate;
1323 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1324 const Scalar control_fraction = fractions[water_pos];
1325 if (control_fraction != 0.0) {
1326 for (
int p = 0; p<np; ++p) {
1327 ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
1333 case Well::ProducerCMode::GRAT:
1335 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1336 Scalar current_rate = -ws.surface_rates[gas_pos];
1339 if (current_rate > 0.0) {
1340 for (
int p = 0; p<np; ++p) {
1341 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1344 const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
1345 const Scalar control_fraction = fractions[gas_pos];
1346 if (control_fraction != 0.0) {
1347 for (
int p = 0; p<np; ++p) {
1348 ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
1356 case Well::ProducerCMode::LRAT:
1358 const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1359 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1360 Scalar current_rate = - ws.surface_rates[water_pos]
1361 - ws.surface_rates[oil_pos];
1364 if (current_rate > 0.0) {
1365 for (
int p = 0; p<np; ++p) {
1366 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1369 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1370 const Scalar control_fraction = fractions[water_pos] + fractions[oil_pos];
1371 if (control_fraction != 0.0) {
1372 for (
int p = 0; p<np; ++p) {
1373 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1379 case Well::ProducerCMode::CRAT:
1382 fmt::format(
"CRAT control not supported, well {}", this->name()),
1385 case Well::ProducerCMode::RESV:
1387 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1388 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1389 Scalar total_res_rate = 0.0;
1390 for (
int p = 0; p<np; ++p) {
1391 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1393 if (controls.prediction_mode) {
1396 if (total_res_rate > 0.0) {
1397 for (
int p = 0; p<np; ++p) {
1398 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1401 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1402 for (
int p = 0; p<np; ++p) {
1403 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1407 std::vector<Scalar> hrates(this->number_of_phases_,0.);
1408 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1409 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1410 hrates[phase_pos] = controls.water_rate;
1412 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1413 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1414 hrates[phase_pos] = controls.oil_rate;
1416 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1417 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1418 hrates[phase_pos] = controls.gas_rate;
1420 std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
1421 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1422 Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1425 if (total_res_rate > 0.0) {
1426 for (
int p = 0; p<np; ++p) {
1427 ws.surface_rates[p] *= target/total_res_rate;
1430 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1431 for (
int p = 0; p<np; ++p) {
1432 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1438 case Well::ProducerCMode::BHP:
1440 ws.bhp = controls.bhp_limit;
1442 for (
int p = 0; p<np; ++p) {
1443 total_rate -= ws.surface_rates[p];
1448 if (total_rate <= 0.0){
1449 for (
int p = 0; p<np; ++p) {
1450 ws.surface_rates[p] = -ws.well_potentials[p];
1455 case Well::ProducerCMode::THP:
1457 const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger);
1459 if (!update_success) {
1463 auto rates = ws.surface_rates;
1464 this->adaptRatesForVFP(rates);
1466 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1468 ws.thp = this->getTHPConstraint(summaryState);
1472 const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1473 if (total_rate <= 0.0) {
1474 for (
int p = 0; p < this->number_of_phases_; ++p) {
1475 ws.surface_rates[p] = -ws.well_potentials[p];
1481 case Well::ProducerCMode::GRUP:
1483 assert(well.isAvailableForGroupControl());
1484 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1485 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1486 well_state[well.name()].efficiency_scaling_factor;
1487 Scalar scale = this->getGroupProductionTargetRate(group,
1497 for (
int p = 0; p<np; ++p) {
1498 ws.surface_rates[p] *= scale;
1500 ws.trivial_group_target =
false;
1504 ws.trivial_group_target =
true;
1508 case Well::ProducerCMode::CMODE_UNDEFINED:
1511 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1517 ws.bhp = controls.bhp_limit;
1522 template<
typename TypeTag>
1531 const bool isGroupControlled = this->wellUnderGroupControl(well_state.
well(this->index_of_well_));
1532 if (!isGroupControlled) {
1534 const auto& summaryState = simulator.vanguard().summaryState();
1535 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1537 return this->wellUnderZeroGroupRateTarget(simulator, well_state, deferred_logger, isGroupControlled);
1541 template <
typename TypeTag>
1546 const std::optional<bool> group_control)
const
1549 const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.
well(this->index_of_well_)));
1550 if (isGroupControlled) {
1551 const auto& summaryState = simulator.vanguard().summaryState();
1552 const auto& group_state = simulator.problem().wellModel().groupState();
1553 const auto& schedule = simulator.vanguard().schedule();
1554 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state, deferred_logger);
1559 template<
typename TypeTag>
1568 return this->wellIsStopped()
1569 || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
1572 template<
typename TypeTag>
1573 std::vector<typename WellInterface<TypeTag>::Scalar>
1579 const int np = this->number_of_phases_;
1580 std::vector<Scalar> scaling_factor(np);
1581 const auto& ws = well_state.
well(this->index_of_well_);
1583 Scalar total_potentials = 0.0;
1584 for (
int p = 0; p<np; ++p) {
1585 total_potentials += ws.well_potentials[p];
1587 if (total_potentials > 0) {
1588 for (
int p = 0; p<np; ++p) {
1589 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1591 return scaling_factor;
1596 const int nperf = this->number_of_local_perforations_;
1597 for (
int perf = 0; perf < nperf; ++perf) {
1598 total_tw += this->well_index_[perf];
1600 total_tw = this->parallelWellInfo().communication().sum(total_tw);
1602 for (
int perf = 0; perf < nperf; ++perf) {
1603 const int cell_idx = this->well_cells_[perf];
1604 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1605 const auto& fs = intQuants.fluidState();
1606 const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
1607 Scalar total_mobility = 0.0;
1608 for (
int p = 0; p < np; ++p) {
1609 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1610 total_mobility += fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value();
1612 for (
int p = 0; p < np; ++p) {
1613 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1614 scaling_factor[p] += well_tw_fraction * fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value() / total_mobility;
1617 return scaling_factor;
1622 template <
typename TypeTag>
1629 assert(this->isProducer());
1633 auto& ws = well_state.
well(this->index_of_well_);
1634 int nonzero_rate_index = -1;
1635 const Scalar floating_point_error_epsilon = 1e-14;
1636 for (
int p = 0; p < this->number_of_phases_; ++p) {
1637 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1638 if (nonzero_rate_index == -1) {
1639 nonzero_rate_index = p;
1648 std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
1649 bool rates_evaluated_at_1bar =
false;
1651 const auto& summary_state = simulator.vanguard().summaryState();
1652 const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
1653 const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
1654 this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
1656 rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
1658 if (std::any_of(well_q_s.begin(), well_q_s.end(), [](
Scalar q) { return q > 0.0; })) {
1660 if (!rates_evaluated_at_1bar) {
1661 this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
1662 rates_evaluated_at_1bar =
true;
1665 for (
auto& q : well_q_s) {
1666 q = std::min(q,
Scalar{0.0});
1671 if (nonzero_rate_index == -1) {
1675 const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
1676 for (
int p = 0; p < this->number_of_phases_; ++p) {
1677 ws.surface_rates[p] = factor * well_q_s[p];
1686 const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1687 const Scalar computed_rate = well_q_s[nonzero_rate_index];
1688 if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {
1690 const Scalar factor = initial_nonzero_rate / computed_rate;
1691 assert(factor < 1.0);
1692 for (
int p = 0; p < this->number_of_phases_; ++p) {
1694 if (p != nonzero_rate_index) {
1695 ws.surface_rates[p] = factor * well_q_s[p];
1704 for (
int p = 0; p < this->number_of_phases_; ++p) {
1705 ws.surface_rates[p] = well_q_s[p];
1709 template <
typename TypeTag>
1710 std::vector<typename WellInterface<TypeTag>::Scalar>
1717 OPM_TIMEFUNCTION_LOCAL();
1720 if (
static_cast<std::size_t
>(perf) >= this->well_cells_.size()) {
1721 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!");
1723 auto wi = std::vector<Scalar>
1724 (this->num_conservation_quantities_, this->well_index_[perf] * trans_mult);
1726 if constexpr (! Indices::gasEnabled) {
1730 const auto& wdfac = this->well_ecl_.getWDFAC();
1732 if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1736 const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
1743 const auto& connection = this->well_ecl_.getConnections()[ws.
perf_data.ecl_index[perf]];
1744 const Scalar Kh = connection.Kh();
1745 const Scalar scaling = 3.141592653589 * Kh * connection.wpimult();
1746 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
1749 const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1750 const Scalar drawdown = cell_pressure - connection_pressure;
1751 const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1752 const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1754 const Scalar b = 2*scaling/wi[gas_comp_idx];
1755 const Scalar c = -2*scaling*mob_g*drawdown;
1757 Scalar consistent_Q = -1.0e20;
1759 const Scalar r2n = b*b + 4*a*c;
1761 const Scalar rn = std::sqrt(r2n);
1762 const Scalar xn1 = (b-rn)*0.5/a;
1766 const Scalar xn2 = (b+rn)*0.5/a;
1767 if (xn2 <= 0 && xn2 > consistent_Q) {
1773 const Scalar r2p = b*b - 4*a*c;
1775 const Scalar rp = std::sqrt(r2p);
1776 const Scalar xp1 = (rp-b)*0.5/a;
1777 if (xp1 > 0 && xp1 < consistent_Q) {
1780 const Scalar xp2 = -(rp+b)*0.5/a;
1781 if (xp2 > 0 && xp2 < consistent_Q) {
1785 wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1790 template <
typename TypeTag>
1796 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1800 auto& d_factor = ws.
perf_data.connection_d_factor;
1802 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1803 const int cell_idx = this->well_cells_[perf];
1804 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1806 d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
1810 template <
typename TypeTag>
1817 auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
1818 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
1822 auto gas_visc = [connection_pressure = ws.
perf_data.pressure[perf],
1824 regIdx = this->pvtRegionIdx(), &intQuants]()
1826 const auto rv = getValue(intQuants.fluidState().Rv());
1828 const auto& gasPvt = FluidSystem::gasPvt();
1833 const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
1834 (regIdx, temperature, connection_pressure);
1836 if (! (rv < rv_sat)) {
1837 return gasPvt.saturatedViscosity(regIdx, temperature,
1838 connection_pressure);
1841 return gasPvt.viscosity(regIdx, temperature, connection_pressure,
1842 rv, getValue(intQuants.fluidState().Rvw()));
1845 const auto& connection = this->well_ecl_.getConnections()
1848 return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
1852 template <
typename TypeTag>
1858 auto connCF = [&connIx = std::as_const(ws.
perf_data.ecl_index),
1859 &conns = this->well_ecl_.getConnections()]
1862 return conns[connIx[perf]].CF();
1865 auto& tmult = ws.
perf_data.connection_compaction_tmult;
1866 auto& ctf = ws.
perf_data.connection_transmissibility_factor;
1868 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1869 const int cell_idx = this->well_cells_[perf];
1871 const auto& intQuants = simulator.model()
1872 .intensiveQuantities(cell_idx, 0);
1874 tmult[perf] = simulator.problem()
1875 .template wellTransMultiplier<double>(intQuants, cell_idx);
1877 ctf[perf] = connCF(perf) * tmult[perf];
1882 template<
typename TypeTag>
1886 if constexpr (Indices::oilEnabled) {
1887 return fs.pressure(FluidSystem::oilPhaseIdx);
1888 }
else if constexpr (Indices::gasEnabled) {
1889 return fs.pressure(FluidSystem::gasPhaseIdx);
1891 return fs.pressure(FluidSystem::waterPhaseIdx);
1895 template <
typename TypeTag>
1896 template<
class Value,
class Callback>
1900 const int local_perf_index,
1901 std::vector<Value>& mob,
1902 Callback& extendEval,
1905 auto relpermArray = []()
1907 if constexpr (std::is_same_v<Value, Scalar>) {
1908 return std::array<Scalar,3>{};
1910 return std::array<Eval,3>{};
1913 if (
static_cast<std::size_t
>(local_perf_index) >= this->well_cells_.size()) {
1914 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!");
1916 const int cell_idx = this->well_cells_[local_perf_index];
1917 assert (
int(mob.size()) == this->num_conservation_quantities_);
1918 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1919 const auto& materialLawManager = simulator.problem().materialLawManager();
1923 const int satid = this->saturation_table_number_[local_perf_index] - 1;
1924 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1925 if (satid == satid_elem) {
1926 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1927 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1931 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
1932 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1934 if constexpr (has_solvent) {
1935 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1938 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1939 auto relativePerms = relpermArray();
1940 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1943 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1946 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1947 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1951 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
1952 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1956 if constexpr (has_solvent) {
1957 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1961 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1962 const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
1963 const auto& connections = this->well_ecl_.getConnections();
1964 const auto& connection = connections[perf_ecl_index];
1965 if (connection.filterCakeActive()) {
1966 std::transform(mob.begin(), mob.end(), mob.begin(),
1967 [mult = this->inj_fc_multiplier_[local_perf_index] ](
const auto val)
1968 { return val * mult; });
1974 template<
typename TypeTag>
1982 const auto& summary_state = simulator.vanguard().summaryState();
1984 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1985 simulator, summary_state, this->getALQ(well_state), deferred_logger,
false);
1986 if (bhp_at_thp_limit) {
1987 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
1988 if (thp_update_iterations) {
1989 computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
1990 rates, deferred_logger);
1992 computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
1993 rates, deferred_logger);
1995 auto& ws = well_state.
well(this->name());
1996 ws.surface_rates = rates;
1997 ws.bhp = *bhp_at_thp_limit;
1998 ws.thp = this->getTHPConstraint(summary_state);
2005 template <
typename TypeTag>
2010 const std::vector<Scalar>& mobility,
2013 const int np = this->number_of_phases_;
2014 for (
int p = 0; p < np; ++p) {
2017 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
2018 const auto connMob =
2019 mobility[FluidSystem::activePhaseToActiveCompIdx(p)] * fs.invB(canonical_phase_idx).value();
2021 connPI[p] = connPICalc(connMob);
2024 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2025 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2027 const auto io = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2028 const auto ig = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2030 const auto vapoil = connPI[ig] * fs.Rv().value();
2031 const auto disgas = connPI[io] * fs.Rs().value();
2033 connPI[io] += vapoil;
2034 connPI[ig] += disgas;
2039 template <
typename TypeTag>
2043 const Phase preferred_phase,
2045 const std::vector<Scalar>& mobility,
2050 if (preferred_phase == Phase::GAS) {
2051 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2053 else if (preferred_phase == Phase::OIL) {
2054 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2056 else if (preferred_phase == Phase::WATER) {
2057 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
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 const int canonicalPhaseIdx = FluidSystem::activeToCanonicalPhaseIdx(phase_pos);
2069 connII[phase_pos] = connIICalc(mt * fs.invB(canonicalPhaseIdx).value());
2072 template<
typename TypeTag>
2073 template<
class GasLiftSingleWell>
2074 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(),
2099 std::set<int> sync_groups;
2100 const auto& summary_state = simulator.vanguard().summaryState();
2101 return std::make_unique<GasLiftSingleWell>(*
this,
#define OPM_DEFLOG_THROW(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:45
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:41
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:198
Definition: SingleWellState.hpp:43
Scalar temperature
Definition: SingleWellState.hpp:99
PerfData< Scalar > perf_data
Definition: SingleWellState.hpp:120
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
Scalar calculateMinimumBhpFromThp(const WellState< Scalar, IndexTraits > &well_state, const Well &well, const SummaryState &summaryState, const Scalar rho) const
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
bool isStableSolution(const WellState< Scalar, IndexTraits > &well_state, const Well &well, const std::vector< Scalar > &rates, const SummaryState &summaryState) const
EvalWell calculateBhpFromThp(const WellState< Scalar, IndexTraits > &well_state, const std::vector< EvalWell > &rates, const Well &well, const SummaryState &summaryState, const Scalar rho, DeferredLogger &deferred_logger) const
std::optional< Scalar > estimateStableBhp(const WellState< Scalar, IndexTraits > &well_state, const Well &well, const std::vector< Scalar > &rates, const Scalar rho, const SummaryState &summaryState) const
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:300
int number_of_local_perforations_
Definition: WellInterfaceGeneric.hpp:336
FluidSystem::Scalar wsolvent_
Definition: WellInterfaceGeneric.hpp:376
const std::vector< FluidSystem::Scalar > & wellIndex() const
Definition: WellInterfaceGeneric.hpp:149
Definition: WellInterfaceIndices.hpp:34
bool solveWellWithBhp(const Simulator &simulator, const double dt, const Scalar bhp, WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:696
void checkWellOperability(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:986
Scalar woxygen() const
Definition: WellInterface_impl.hpp:165
IndividualOrGroup
Definition: WellInterface.hpp:244
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:81
Scalar computeConnectionDFactor(const int perf, const IntensiveQuantities &intQuants, const SingleWellStateType &ws) const
Definition: WellInterface_impl.hpp:1813
void solveWellEquation(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:776
bool updateWellControlAndStatusLocalIteration(const Simulator &simulator, WellStateType &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:272
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:103
bool iterateWellEquations(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:500
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:135
Scalar wfoam() const
Definition: WellInterface_impl.hpp:127
std::optional< Scalar > estimateOperableBhp(const Simulator &ebos_simulator, const double dt, WellStateType &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:663
std::vector< RateVector > connectionRates_
Definition: WellInterface.hpp:367
void computeConnLevelProdInd(const FluidState &fs, const std::function< Scalar(const Scalar)> &connPICalc, const std::vector< Scalar > &mobility, Scalar *connPI) const
Definition: WellInterface_impl.hpp:2008
void addCellRates(RateVector &rates, int cellIdx) const
Definition: WellInterface_impl.hpp:950
void gliftBeginTimeStepWellTestUpdateALQ(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, GLiftEclWells &ecl_well_map, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1013
void assembleWellEq(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:838
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const
Definition: WellInterface_impl.hpp:966
virtual void init(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 getMobility(const Simulator &simulator, const int local_perf_index, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1899
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:854
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:86
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:82
std::vector< Scalar > initialWellRateFractions(const Simulator &ebosSimulator, const WellStateType &well_state) const
Definition: WellInterface_impl.hpp:1575
bool updateWellStateWithTHPTargetProd(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1977
void updateConnectionDFactor(const Simulator &simulator, SingleWellStateType &ws) const
Definition: WellInterface_impl.hpp:1793
Eval getPerfCellPressure(const FluidState &fs) const
Definition: WellInterface_impl.hpp:1884
void initializeProducerWellState(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1625
bool wellUnderZeroGroupRateTarget(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger, std::optional< bool > group_control=std::nullopt) const
Definition: WellInterface_impl.hpp:1543
typename Base::ModelParameters ModelParameters
Definition: WellInterface.hpp:109
bool updateWellControl(const Simulator &simulator, const IndividualOrGroup iog, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:189
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:83
static constexpr bool has_solvent
Definition: WellInterface.hpp:111
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: WellInterface.hpp:89
void updateConnectionTransmissibilityFactor(const Simulator &simulator, SingleWellStateType &ws) const
Definition: WellInterface_impl.hpp:1855
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:2042
typename GasLiftGroupInfo< Scalar, IndexTraits >::GLiftEclWells GLiftEclWells
Definition: WellInterface.hpp:91
std::unique_ptr< GasLiftSingleWell > initializeGliftWellTest_(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, GLiftEclWells &ecl_well_map, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:2076
void wellTesting(const Simulator &simulator, const double simulation_time, WellStateType &well_state, const GroupState< Scalar > &group_state, WellTestState &welltest_state, GLiftEclWells &ecl_well_map, std::map< std::string, double > &open_times, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:380
Scalar wsalt() const
Definition: WellInterface_impl.hpp:141
void prepareWellBeforeAssembling(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:874
typename Base::Eval Eval
Definition: WellInterface.hpp:95
bool wellUnderZeroRateTarget(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1525
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_conservation_quantities, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Constructor.
Definition: WellInterface_impl.hpp:58
bool solveWellForTesting(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:754
void updateWellOperability(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1075
bool solveWellWithOperabilityCheck(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:563
Scalar wpolymer() const
Definition: WellInterface_impl.hpp:111
bool updateWellOperabilityFromWellEq(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1112
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:85
Scalar wurea() const
Definition: WellInterface_impl.hpp:177
bool solveWellWithZeroRate(const Simulator &simulator, const double dt, WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:733
Scalar wmicrobes() const
Definition: WellInterface_impl.hpp:153
virtual void updateWellStateWithTarget(const Simulator &simulator, const GroupState< Scalar > &group_state, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1131
static constexpr bool has_zFraction
Definition: WellInterface.hpp:112
bool stoppedOrZeroRateTarget(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1562
Definition: WellState.hpp:66
constexpr int numPhases() const
The number of phases present.
Definition: WellState.hpp:248
const SingleWellState< Scalar, IndexTraits > & well(std::size_t well_index) const
Definition: WellState.hpp:289
@ NONE
Definition: DeferredLogger.hpp:46
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:39
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Static data associated with a well perforation.
Definition: PerforationData.hpp:30