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 bool is_grup =
false;
206 if (well.isInjector()) {
207 from = WellInjectorCMode2String(ws.injection_cmode);
208 is_grup = ws.injection_cmode == Well::InjectorCMode::GRUP;
210 from = WellProducerCMode2String(ws.production_cmode);
211 is_grup = ws.production_cmode == Well::ProducerCMode::GRUP;
214 const int episodeIdx = simulator.episodeIndex();
215 const int iterationIdx = simulator.model().newtonMethod().numIterations();
216 const int nupcol = schedule[episodeIdx].nupcol();
217 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
218 if (oscillating && !is_grup) {
220 const bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
222 const auto msg = fmt::format(
" The control mode for well {} is oscillating. \n"
223 "We don't allow for more than {} switches after NUPCOL iterations. (NUPCOL = {}) \n"
224 "The control is kept at {}.",
225 this->name(), this->param_.max_number_of_well_switches_, nupcol, from);
226 deferred_logger.
info(msg);
228 this->well_control_log_.push_back(from);
232 bool changed =
false;
233 if (iog == IndividualOrGroup::Individual) {
234 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
235 }
else if (iog == IndividualOrGroup::Group) {
236 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState,
true, deferred_logger);
238 assert(iog == IndividualOrGroup::Both);
239 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
245 if (well.isInjector()) {
246 to = WellInjectorCMode2String(ws.injection_cmode);
248 to = WellProducerCMode2String(ws.production_cmode);
250 std::ostringstream ss;
251 ss <<
" Switching control mode for well " << this->name()
255 ss <<
" on rank " << cc.rank();
257 deferred_logger.
debug(ss.str());
262 if (iterationIdx >= nupcol || this->well_control_log_.empty()) {
263 this->well_control_log_.push_back(from);
265 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
266 updatePrimaryVariables(simulator, well_state, deferred_logger);
272 template<
typename TypeTag>
278 const Well::InjectionControls& inj_controls,
279 const Well::ProductionControls& prod_controls,
282 const bool fixed_control,
283 const bool fixed_status)
286 const auto& summary_state = simulator.vanguard().summaryState();
287 const auto& schedule = simulator.vanguard().schedule();
288 auto& ws = well_state.
well(this->index_of_well_);
290 if (this->isInjector()) {
291 from = WellInjectorCMode2String(ws.injection_cmode);
293 from = WellProducerCMode2String(ws.production_cmode);
295 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
297 if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(well_state.
well(this->index_of_well_).status == WellStatus::OPEN)) {
301 const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
302 if (!this->wellIsStopped()){
303 if (wqTotal*sgn <= 0.0 && !fixed_status){
307 bool changed =
false;
308 if (!fixed_control) {
313 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
314 prod_controls.hasControl(Well::ProducerCMode::GRUP);
315 bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
316 if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
317 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
319 if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
320 changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,
false, deferred_logger);
324 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
325 ws.production_cmode == Well::ProducerCMode::THP;
327 ws.thp = this->getTHPConstraint(summary_state);
330 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
332 updatePrimaryVariables(simulator, well_state, deferred_logger);
337 }
else if (!fixed_status){
339 const Scalar bhp = well_state.
well(this->index_of_well_).bhp;
340 Scalar prod_limit = prod_controls.bhp_limit;
341 Scalar inj_limit = inj_controls.bhp_limit;
342 const bool has_thp = this->wellHasTHPConstraints(summary_state);
344 std::vector<Scalar> rates(this->num_conservation_quantities_);
345 if (this->isInjector()){
347 calculateBhpFromThp(well_state, rates,
350 this->getRefDensity(),
352 inj_limit = std::min(bhp_thp,
static_cast<Scalar>(inj_controls.bhp_limit));
357 calculateMinimumBhpFromThp(well_state,
360 this->getRefDensity());
361 prod_limit = std::max(bhp_min,
static_cast<Scalar>(prod_controls.bhp_limit));
364 const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
367 well_state.
well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
369 well_state.
well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
380 template<
typename TypeTag>
384 const double simulation_time,
387 WellTestState& well_test_state,
389 std::map<std::string, double>& open_times,
393 deferred_logger.
info(
" well " + this->name() +
" is being tested");
396 auto& ws = well_state_copy.
well(this->indexOfWell());
398 const auto& summary_state = simulator.vanguard().summaryState();
399 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
400 if (this->isProducer()) {
401 ws.production_cmode = has_thp_limit ? Well::ProducerCMode::THP : Well::ProducerCMode::BHP;
403 ws.injection_cmode = has_thp_limit ? Well::InjectorCMode::THP : Well::InjectorCMode::BHP;
408 scaleSegmentRatesAndPressure(well_state_copy);
409 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
410 updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
412 if (this->isProducer()) {
413 const auto& schedule = simulator.vanguard().schedule();
414 const auto report_step = simulator.episodeIndex();
415 const auto& glo = schedule.glo(report_step);
417 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() && (well_state.
well(this->index_of_well_).status == 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 converged =
true;
577 auto& ws = well_state.
well(this->index_of_well_);
579 if (this->wellIsStopped()) {
581 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
582 if (!bhp_target.has_value()) {
584 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
585 deferred_logger.
debug(msg);
588 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
590 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
591 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
595 ws.thp = this->getTHPConstraint(summary_state);
596 const Scalar bhp = std::max(bhp_target.value(),
597 static_cast<Scalar>(prod_controls.bhp_limit));
598 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
602 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
605 const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
607 if (converged && !stoppedOrZeroRateTarget(simulator, well_state, deferred_logger) && isThp) {
608 auto rates = well_state.
well(this->index_of_well_).surface_rates;
609 this->adaptRatesForVFP(rates);
610 this->updateIPRImplicit(simulator, well_state, deferred_logger);
614 this->operability_status_.use_vfpexplicit =
true;
617 const Scalar reltol = 1e-3;
618 const Scalar cur_bhp = ws.bhp;
619 if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
620 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
621 deferred_logger.
debug(msg);
622 solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger);
624 ws.thp = this->getTHPConstraint(summary_state);
625 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
632 this->operability_status_.use_vfpexplicit =
true;
634 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
635 if (!bhp_target.has_value()) {
638 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
640 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
641 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
645 const Scalar bhp = std::max(bhp_target.value(),
646 static_cast<Scalar>(prod_controls.bhp_limit));
647 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
648 ws.thp = this->getTHPConstraint(summary_state);
649 const auto msg = fmt::format(
"Well {} did not converge, re-solving with explicit fractions for VFP caculations.", this->name());
650 deferred_logger.
debug(msg);
651 converged = this->iterateWellEqWithSwitching(simulator, dt,
660 this->operability_status_.can_obtain_bhp_with_thp_limit = !this->wellIsStopped();
661 this->operability_status_.obey_thp_limit_under_bhp_limit = !this->wellIsStopped();
665 template<
typename TypeTag>
666 std::optional<typename WellInterface<TypeTag>::Scalar>
671 const SummaryState& summary_state,
674 if (!this->wellHasTHPConstraints(summary_state)) {
676 const bool converged = solveWellWithBhp(simulator, dt, bhp_limit, well_state, deferred_logger);
677 if (!converged || this->wellIsStopped()) {
688 const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
689 if (!converged || this->wellIsStopped()) {
692 this->updateIPRImplicit(simulator, well_state, deferred_logger);
693 auto rates = well_state.
well(this->index_of_well_).surface_rates;
694 this->adaptRatesForVFP(rates);
698 template<
typename TypeTag>
710 auto inj_controls = Well::InjectionControls(0);
711 auto prod_controls = Well::ProductionControls(0);
712 auto& ws = well_state.
well(this->index_of_well_);
713 auto cmode_inj = ws.injection_cmode;
714 auto cmode_prod = ws.production_cmode;
715 if (this->isInjector()) {
716 inj_controls.addControl(Well::InjectorCMode::BHP);
717 inj_controls.bhp_limit = bhp;
718 inj_controls.cmode = Well::InjectorCMode::BHP;
719 ws.injection_cmode = Well::InjectorCMode::BHP;
721 prod_controls.addControl(Well::ProducerCMode::BHP);
722 prod_controls.bhp_limit = bhp;
723 prod_controls.cmode = Well::ProducerCMode::BHP;
724 ws.production_cmode = Well::ProducerCMode::BHP;
729 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true);
730 ws.injection_cmode = cmode_inj;
731 ws.production_cmode = cmode_prod;
735 template<
typename TypeTag>
745 const auto well_status_orig = this->wellStatus_;
749 auto inj_controls = Well::InjectionControls(0);
750 auto prod_controls = Well::ProductionControls(0);
751 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true,
true);
752 this->wellStatus_ = well_status_orig;
756 template<
typename TypeTag>
765 const double dt = simulator.timeStepSize();
767 const auto& summary_state = simulator.vanguard().summaryState();
768 auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
769 auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
770 this->onlyKeepBHPandTHPcontrols(summary_state, well_state, inj_controls, prod_controls);
772 bool converged =
false;
775 if (!this->param_.local_well_solver_control_switching_){
776 converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
778 if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.
well(this->index_of_well_).status == WellStatus::OPEN)) {
779 converged = solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
781 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
785 }
catch (NumericalProblem& e ) {
786 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
787 deferred_logger.
warning(
"INNER_ITERATION_FAILED", msg);
792 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" converged");
795 const int max_iter = this->param_.max_welleq_iter_;
796 deferred_logger.
debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
802 template<
typename TypeTag>
811 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
816 const double dt = simulator.timeStepSize();
817 bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
827 auto& ws = well_state.
well(this->indexOfWell());
828 bool thp_control =
false;
829 if (this->well_ecl_.isInjector()) {
830 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
832 ws.injection_cmode = Well::InjectorCMode::BHP;
833 if (this->well_control_log_.empty()) {
834 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
838 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
840 ws.production_cmode = Well::ProducerCMode::BHP;
841 if (this->well_control_log_.empty()) {
842 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
847 const std::string msg = std::string(
"The newly opened well ") + this->name()
848 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
849 deferred_logger.
debug(msg);
850 converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
855 const int max_iter = this->param_.max_welleq_iter_;
856 deferred_logger.
debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
858 well_state = well_state0;
864 template <
typename TypeTag>
874 prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
875 assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
880 template <
typename TypeTag>
890 const auto& summary_state = simulator.vanguard().summaryState();
891 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
892 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
895 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
900 template<
typename TypeTag>
910 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
912 if (this->param_.check_well_operability_iter_)
913 checkWellOperability(simulator, well_state, deferred_logger);
916 const int iteration_idx = simulator.model().newtonMethod().numIterations();
917 if (iteration_idx < this->param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
918 const auto& ws = well_state.
well(this->indexOfWell());
919 const bool nonzero_rate_original =
920 std::any_of(ws.surface_rates.begin(),
921 ws.surface_rates.begin() + well_state.
numPhases(),
922 [](
Scalar rate) { return rate != Scalar(0.0); });
924 this->operability_status_.solvable =
true;
925 if (number_of_well_reopenings_ >= this->param_.max_well_status_switch_) {
927 if (number_of_well_reopenings_ == this->param_.max_well_status_switch_) {
928 const std::string msg = fmt::format(
"well {} is oscillating between open and stop. \n"
929 "We don't allow for more than {} re-openings "
930 "and the well is therefore kept stopped.",
931 this->name(), number_of_well_reopenings_);
932 deferred_logger.
debug(msg);
935 changed_to_stopped_this_step_ =
true;
936 bool converged_zero_rate = this->solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
937 if (this->param_.shut_unsolvable_wells_ && !converged_zero_rate ) {
938 this->operability_status_.solvable =
false;
941 number_of_well_reopenings_++;
944 bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
947 const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
948 if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
952 this->operability_status_.resetOperability();
954 deferred_logger.
debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
955 number_of_well_reopenings_++;
959 if (this->param_.shut_unsolvable_wells_) {
960 this->operability_status_.solvable =
false;
964 if (this->operability_status_.has_negative_potentials) {
965 auto well_state_copy = well_state;
966 std::vector<Scalar> potentials;
968 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
969 }
catch (
const std::exception& e) {
970 const std::string msg = fmt::format(
"well {}: computeWellPotentials() failed "
971 "during attempt to recompute potentials for well: ",
972 this->name(), e.what());
973 deferred_logger.
info(msg);
974 this->operability_status_.has_negative_potentials =
true;
976 auto& ws = well_state.
well(this->indexOfWell());
978 for (
int p = 0; p < np; ++p) {
979 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
982 this->changed_to_open_this_step_ =
false;
983 changed_to_stopped_this_step_ =
false;
985 const bool well_operable = this->operability_status_.isOperableAndSolvable();
986 if (!well_operable) {
988 if (old_well_operable) {
989 deferred_logger.
debug(
" well " + this->name() +
" gets STOPPED during iteration ");
990 changed_to_stopped_this_step_ =
true;
992 }
else if (well_state.
isOpen(this->name())) {
994 if (!old_well_operable) {
995 deferred_logger.
debug(
" well " + this->name() +
" gets REVIVED during iteration ");
996 this->changed_to_open_this_step_ =
true;
1001 template<
typename TypeTag>
1005 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
1008 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
1009 if (this->cells()[perfIdx] == cellIdx) {
1010 for (
int i = 0; i < RateVector::dimension; ++i) {
1011 rates[i] += connectionRates_[perfIdx][i];
1017 template<
typename TypeTag>
1021 for (
int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
1022 if (this->cells()[perfIdx] == cellIdx) {
1023 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
1024 return connectionRates_[perfIdx][activeCompIdx].value();
1028 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
1036 template<
typename TypeTag>
1044 if (!this->param_.check_well_operability_) {
1048 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
1052 updateWellOperability(simulator, well_state, deferred_logger);
1053 if (!this->operability_status_.isOperableAndSolvable()) {
1054 this->operability_status_.use_vfpexplicit =
true;
1055 deferred_logger.
debug(
"EXPLICIT_LOOKUP_VFP",
1056 "well not operable, trying with explicit vfp lookup: " + this->name());
1057 updateWellOperability(simulator, well_state, deferred_logger);
1063 template<
typename TypeTag>
1073 const auto& summary_state = simulator.vanguard().summaryState();
1074 const auto& well_name = this->name();
1075 if (!this->wellHasTHPConstraints(summary_state)) {
1076 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1077 deferred_logger.
info(msg);
1080 const auto& schedule = simulator.vanguard().schedule();
1081 const auto report_step_idx = simulator.episodeIndex();
1082 const auto& glo = schedule.glo(report_step_idx);
1083 if (!glo.has_well(well_name)) {
1084 const std::string msg = fmt::format(
1085 "GLIFT WTEST: Well {} : Gas lift not activated: "
1086 "WLIFTOPT is probably missing. Skipping.", well_name);
1087 deferred_logger.
info(msg);
1090 const auto& gl_well = glo.well(well_name);
1093 std::unique_ptr<GasLiftSingleWell> glift =
1094 initializeGliftWellTest_<GasLiftSingleWell>(simulator,
1099 auto [wtest_alq, success] = glift->wellTestALQ();
1101 const auto& unit_system = schedule.getUnits();
1103 well_state.
well(well_name).alq_state.set(wtest_alq);
1105 "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1106 well_name, unit_system.from_si(UnitSystem::measure::gas_surface_rate, wtest_alq));
1109 if (!gl_well.use_glo()) {
1111 "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1113 unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.
well(well_name).alq_state.get()));
1118 "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1122 deferred_logger.
info(msg);
1125 template<
typename TypeTag>
1133 if (this->param_.local_well_solver_control_switching_) {
1134 const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger);
1136 this->operability_status_.solvable =
false;
1137 deferred_logger.
debug(
"Operability check using well equations did not converge for well "
1138 + this->name() +
". Mark the well as unsolvable." );
1142 this->operability_status_.resetOperability();
1144 bool thp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1145 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1146 bool bhp_controlled = this->isInjector() ? well_state.
well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1147 well_state.
well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1151 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
1152 if (check_thp || bhp_controlled) {
1153 updateIPR(simulator, deferred_logger);
1154 checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
1158 checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger);
1162 template<
typename TypeTag>
1171 assert(this->param_.local_well_solver_control_switching_);
1172 this->operability_status_.resetOperability();
1174 const auto& group_state = simulator.problem().wellModel().groupState();
1175 const double dt = simulator.timeStepSize();
1177 bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger);
1181 template<
typename TypeTag>
1189 template<
typename TypeTag>
1199 const auto& well = this->well_ecl_;
1200 const int well_index = this->index_of_well_;
1201 auto& ws = well_state.
well(well_index);
1203 const auto& summaryState = simulator.vanguard().summaryState();
1204 const auto& schedule = simulator.vanguard().schedule();
1208 ws.primaryvar.resize(0);
1210 if (this->wellIsStopped()) {
1211 for (
int p = 0; p<np; ++p) {
1212 ws.surface_rates[p] = 0;
1218 if (this->isInjector() )
1220 const auto& controls = well.injectionControls(summaryState);
1222 InjectorType injectorType = controls.injector_type;
1224 switch (injectorType) {
1225 case InjectorType::WATER:
1227 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1230 case InjectorType::OIL:
1232 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1235 case InjectorType::GAS:
1237 phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1241 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1244 const auto current = ws.injection_cmode;
1247 case Well::InjectorCMode::RATE:
1249 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1250 if(this->rsRvInj() > 0) {
1251 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1252 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1253 ws.surface_rates[gas_pos] = controls.surface_rate * this->rsRvInj();
1254 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1255 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1256 ws.surface_rates[oil_pos] = controls.surface_rate * this->rsRvInj();
1258 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 );
1264 case Well::InjectorCMode::RESV:
1266 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1267 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
1268 const Scalar coeff = convert_coeff[phasePos];
1269 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1273 case Well::InjectorCMode::THP:
1275 auto rates = ws.surface_rates;
1280 this->getRefDensity(),
1283 ws.thp = this->getTHPConstraint(summaryState);
1288 Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1289 if (total_rate <= 0.0)
1290 ws.surface_rates = ws.well_potentials;
1294 case Well::InjectorCMode::BHP:
1296 ws.bhp = controls.bhp_limit;
1298 for (
int p = 0; p<np; ++p) {
1299 total_rate += ws.surface_rates[p];
1304 if (total_rate <= 0.0)
1305 ws.surface_rates = ws.well_potentials;
1309 case Well::InjectorCMode::GRUP:
1311 assert(well.isAvailableForGroupControl());
1312 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1313 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1314 well_state[well.name()].efficiency_scaling_factor;
1315 std::optional<Scalar> target =
1316 this->getGroupInjectionTargetRate(group,
1325 ws.surface_rates[phasePos] = *target;
1328 case Well::InjectorCMode::CMODE_UNDEFINED:
1330 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
1340 ws.surface_rates[phasePos] = std::max(
Scalar{1.e-7}, ws.surface_rates[phasePos]);
1343 ws.bhp = controls.bhp_limit;
1349 const auto current = ws.production_cmode;
1350 const auto& controls = well.productionControls(summaryState);
1352 case Well::ProducerCMode::ORAT:
1354 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1355 Scalar current_rate = -ws.surface_rates[oil_pos];
1358 if (current_rate > 0.0) {
1359 for (
int p = 0; p<np; ++p) {
1360 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1363 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1364 double control_fraction = fractions[oil_pos];
1365 if (control_fraction != 0.0) {
1366 for (
int p = 0; p<np; ++p) {
1367 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1373 case Well::ProducerCMode::WRAT:
1375 const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1376 Scalar current_rate = -ws.surface_rates[water_pos];
1379 if (current_rate > 0.0) {
1380 for (
int p = 0; p<np; ++p) {
1381 ws.surface_rates[p] *= controls.water_rate/current_rate;
1384 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1385 const Scalar control_fraction = fractions[water_pos];
1386 if (control_fraction != 0.0) {
1387 for (
int p = 0; p<np; ++p) {
1388 ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
1394 case Well::ProducerCMode::GRAT:
1396 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1397 Scalar current_rate = -ws.surface_rates[gas_pos];
1400 if (current_rate > 0.0) {
1401 for (
int p = 0; p<np; ++p) {
1402 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1405 const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
1406 const Scalar control_fraction = fractions[gas_pos];
1407 if (control_fraction != 0.0) {
1408 for (
int p = 0; p<np; ++p) {
1409 ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
1417 case Well::ProducerCMode::LRAT:
1419 const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1420 const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1421 Scalar current_rate = - ws.surface_rates[water_pos]
1422 - ws.surface_rates[oil_pos];
1425 if (current_rate > 0.0) {
1426 for (
int p = 0; p<np; ++p) {
1427 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1430 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1431 const Scalar control_fraction = fractions[water_pos] + fractions[oil_pos];
1432 if (control_fraction != 0.0) {
1433 for (
int p = 0; p<np; ++p) {
1434 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1440 case Well::ProducerCMode::CRAT:
1443 fmt::format(
"CRAT control not supported, well {}", this->name()),
1446 case Well::ProducerCMode::RESV:
1448 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1449 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1450 Scalar total_res_rate = 0.0;
1451 for (
int p = 0; p<np; ++p) {
1452 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1454 if (controls.prediction_mode) {
1457 if (total_res_rate > 0.0) {
1458 for (
int p = 0; p<np; ++p) {
1459 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1462 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1463 for (
int p = 0; p<np; ++p) {
1464 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1468 std::vector<Scalar> hrates(this->number_of_phases_,0.);
1469 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1470 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1471 hrates[phase_pos] = controls.water_rate;
1473 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1474 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1475 hrates[phase_pos] = controls.oil_rate;
1477 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1478 const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1479 hrates[phase_pos] = controls.gas_rate;
1481 std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
1482 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1483 Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1486 if (total_res_rate > 0.0) {
1487 for (
int p = 0; p<np; ++p) {
1488 ws.surface_rates[p] *= target/total_res_rate;
1491 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1492 for (
int p = 0; p<np; ++p) {
1493 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1499 case Well::ProducerCMode::BHP:
1501 ws.bhp = controls.bhp_limit;
1503 for (
int p = 0; p<np; ++p) {
1504 total_rate -= ws.surface_rates[p];
1509 if (total_rate <= 0.0){
1510 for (
int p = 0; p<np; ++p) {
1511 ws.surface_rates[p] = -ws.well_potentials[p];
1516 case Well::ProducerCMode::THP:
1518 const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger);
1520 if (!update_success) {
1524 auto rates = ws.surface_rates;
1525 this->adaptRatesForVFP(rates);
1527 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1529 ws.thp = this->getTHPConstraint(summaryState);
1533 const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1534 if (total_rate <= 0.0) {
1535 for (
int p = 0; p < this->number_of_phases_; ++p) {
1536 ws.surface_rates[p] = -ws.well_potentials[p];
1542 case Well::ProducerCMode::GRUP:
1544 assert(well.isAvailableForGroupControl());
1545 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1546 const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1547 well_state[well.name()].efficiency_scaling_factor;
1548 Scalar scale = this->getGroupProductionTargetRate(group,
1558 for (
int p = 0; p<np; ++p) {
1559 ws.surface_rates[p] *= scale;
1561 ws.trivial_group_target =
false;
1565 ws.trivial_group_target =
true;
1569 case Well::ProducerCMode::CMODE_UNDEFINED:
1572 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1578 ws.bhp = controls.bhp_limit;
1583 template<
typename TypeTag>
1592 const bool isGroupControlled = this->wellUnderGroupControl(well_state.
well(this->index_of_well_));
1593 if (!isGroupControlled) {
1595 const auto& summaryState = simulator.vanguard().summaryState();
1596 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1598 return this->wellUnderZeroGroupRateTarget(simulator, well_state, deferred_logger, isGroupControlled);
1602 template <
typename TypeTag>
1607 const std::optional<bool> group_control)
const
1610 const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.
well(this->index_of_well_)));
1611 if (isGroupControlled) {
1612 const auto& summaryState = simulator.vanguard().summaryState();
1613 const auto& group_state = simulator.problem().wellModel().groupState();
1614 const auto& schedule = simulator.vanguard().schedule();
1615 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state, deferred_logger);
1620 template<
typename TypeTag>
1629 return this->wellIsStopped()
1630 || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
1633 template<
typename TypeTag>
1634 std::vector<typename WellInterface<TypeTag>::Scalar>
1640 const int np = this->number_of_phases_;
1641 std::vector<Scalar> scaling_factor(np);
1642 const auto& ws = well_state.
well(this->index_of_well_);
1644 Scalar total_potentials = 0.0;
1645 for (
int p = 0; p<np; ++p) {
1646 total_potentials += ws.well_potentials[p];
1648 if (total_potentials > 0) {
1649 for (
int p = 0; p<np; ++p) {
1650 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1652 return scaling_factor;
1657 const int nperf = this->number_of_local_perforations_;
1658 for (
int perf = 0; perf < nperf; ++perf) {
1659 total_tw += this->well_index_[perf];
1661 total_tw = this->parallelWellInfo().communication().sum(total_tw);
1663 for (
int perf = 0; perf < nperf; ++perf) {
1664 const int cell_idx = this->well_cells_[perf];
1665 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1666 const auto& fs = intQuants.fluidState();
1667 const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
1668 Scalar total_mobility = 0.0;
1669 for (
int p = 0; p < np; ++p) {
1670 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1671 total_mobility += fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value();
1673 for (
int p = 0; p < np; ++p) {
1674 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1675 scaling_factor[p] += well_tw_fraction * fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value() / total_mobility;
1678 return scaling_factor;
1683 template <
typename TypeTag>
1690 assert(this->isProducer());
1694 auto& ws = well_state.
well(this->index_of_well_);
1695 int nonzero_rate_index = -1;
1696 const Scalar floating_point_error_epsilon = 1e-14;
1697 for (
int p = 0; p < this->number_of_phases_; ++p) {
1698 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1699 if (nonzero_rate_index == -1) {
1700 nonzero_rate_index = p;
1709 std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
1710 bool rates_evaluated_at_1bar =
false;
1712 const auto& summary_state = simulator.vanguard().summaryState();
1713 const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
1714 const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
1715 this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
1717 rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
1719 if (std::any_of(well_q_s.begin(), well_q_s.end(), [](
Scalar q) { return q > 0.0; })) {
1721 if (!rates_evaluated_at_1bar) {
1722 this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
1723 rates_evaluated_at_1bar =
true;
1726 for (
auto& q : well_q_s) {
1727 q = std::min(q,
Scalar{0.0});
1732 if (nonzero_rate_index == -1) {
1736 const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
1737 for (
int p = 0; p < this->number_of_phases_; ++p) {
1738 ws.surface_rates[p] = factor * well_q_s[p];
1747 const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1748 const Scalar computed_rate = well_q_s[nonzero_rate_index];
1749 if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {
1751 const Scalar factor = initial_nonzero_rate / computed_rate;
1752 assert(factor < 1.0);
1753 for (
int p = 0; p < this->number_of_phases_; ++p) {
1755 if (p != nonzero_rate_index) {
1756 ws.surface_rates[p] = factor * well_q_s[p];
1765 for (
int p = 0; p < this->number_of_phases_; ++p) {
1766 ws.surface_rates[p] = well_q_s[p];
1770 template <
typename TypeTag>
1771 std::vector<typename WellInterface<TypeTag>::Scalar>
1778 OPM_TIMEFUNCTION_LOCAL(Subsystem::Wells);
1781 if (
static_cast<std::size_t
>(perf) >= this->well_cells_.size()) {
1782 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!");
1784 auto wi = std::vector<Scalar>
1785 (this->num_conservation_quantities_, this->well_index_[perf] * trans_mult);
1787 if constexpr (! Indices::gasEnabled) {
1791 const auto& wdfac = this->well_ecl_.getWDFAC();
1793 if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1797 const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
1804 const auto& connection = this->well_ecl_.getConnections()[ws.
perf_data.ecl_index[perf]];
1805 const Scalar Kh = connection.Kh();
1806 const Scalar scaling = 3.141592653589 * Kh * connection.wpimult();
1807 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
1810 const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1811 const Scalar drawdown = cell_pressure - connection_pressure;
1812 const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1813 const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1815 const Scalar b = 2*scaling/wi[gas_comp_idx];
1816 const Scalar c = -2*scaling*mob_g*drawdown;
1818 Scalar consistent_Q = -1.0e20;
1820 const Scalar r2n = b*b + 4*a*c;
1822 const Scalar rn = std::sqrt(r2n);
1823 const Scalar xn1 = (b-rn)*0.5/a;
1827 const Scalar xn2 = (b+rn)*0.5/a;
1828 if (xn2 <= 0 && xn2 > consistent_Q) {
1834 const Scalar r2p = b*b - 4*a*c;
1836 const Scalar rp = std::sqrt(r2p);
1837 const Scalar xp1 = (rp-b)*0.5/a;
1838 if (xp1 > 0 && xp1 < consistent_Q) {
1841 const Scalar xp2 = -(rp+b)*0.5/a;
1842 if (xp2 > 0 && xp2 < consistent_Q) {
1846 wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1851 template <
typename TypeTag>
1857 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1861 auto& d_factor = ws.
perf_data.connection_d_factor;
1863 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1864 const int cell_idx = this->well_cells_[perf];
1865 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1867 d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
1871 template <
typename TypeTag>
1878 auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
1879 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
1883 auto gas_visc = [connection_pressure = ws.
perf_data.pressure[perf],
1885 regIdx = this->pvtRegionIdx(), &intQuants]()
1887 const auto rv = getValue(intQuants.fluidState().Rv());
1889 const auto& gasPvt = FluidSystem::gasPvt();
1894 const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
1895 (regIdx, temperature, connection_pressure);
1897 if (! (rv < rv_sat)) {
1898 return gasPvt.saturatedViscosity(regIdx, temperature,
1899 connection_pressure);
1902 return gasPvt.viscosity(regIdx, temperature, connection_pressure,
1903 rv, getValue(intQuants.fluidState().Rvw()));
1906 const auto& connection = this->well_ecl_.getConnections()
1909 return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
1913 template <
typename TypeTag>
1919 auto connCF = [&connIx = std::as_const(ws.
perf_data.ecl_index),
1920 &conns = this->well_ecl_.getConnections()]
1923 return conns[connIx[perf]].CF();
1926 auto& tmult = ws.
perf_data.connection_compaction_tmult;
1927 auto& ctf = ws.
perf_data.connection_transmissibility_factor;
1929 for (
int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1930 const int cell_idx = this->well_cells_[perf];
1932 const auto& intQuants = simulator.model()
1933 .intensiveQuantities(cell_idx, 0);
1935 tmult[perf] = simulator.problem()
1936 .template wellTransMultiplier<double>(intQuants, cell_idx);
1938 ctf[perf] = connCF(perf) * tmult[perf];
1943 template<
typename TypeTag>
1947 if constexpr (Indices::oilEnabled) {
1948 return fs.pressure(FluidSystem::oilPhaseIdx);
1949 }
else if constexpr (Indices::gasEnabled) {
1950 return fs.pressure(FluidSystem::gasPhaseIdx);
1952 return fs.pressure(FluidSystem::waterPhaseIdx);
1956 template <
typename TypeTag>
1957 template<
class Value,
class Callback>
1961 const int local_perf_index,
1962 std::vector<Value>& mob,
1963 Callback& extendEval,
1966 auto relpermArray = []()
1968 if constexpr (std::is_same_v<Value, Scalar>) {
1969 return std::array<Scalar,3>{};
1971 return std::array<Eval,3>{};
1974 if (
static_cast<std::size_t
>(local_perf_index) >= this->well_cells_.size()) {
1975 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!");
1977 const int cell_idx = this->well_cells_[local_perf_index];
1978 assert (
int(mob.size()) == this->num_conservation_quantities_);
1979 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1980 const auto& materialLawManager = simulator.problem().materialLawManager();
1984 const int satid = this->saturation_table_number_[local_perf_index] - 1;
1985 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1986 if (satid == satid_elem) {
1987 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1988 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1992 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
1993 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1995 if constexpr (has_solvent) {
1996 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1999 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
2000 auto relativePerms = relpermArray();
2001 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
2004 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
2007 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2008 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2012 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
2013 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
2017 if constexpr (has_solvent) {
2018 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
2022 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
2023 const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
2024 const auto& connections = this->well_ecl_.getConnections();
2025 const auto& connection = connections[perf_ecl_index];
2026 if (connection.filterCakeActive()) {
2027 std::transform(mob.begin(), mob.end(), mob.begin(),
2028 [mult = this->inj_fc_multiplier_[local_perf_index] ](
const auto val)
2029 { return val * mult; });
2035 template<
typename TypeTag>
2043 const auto& summary_state = simulator.vanguard().summaryState();
2045 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
2046 simulator, summary_state, this->getALQ(well_state), deferred_logger,
false);
2047 if (bhp_at_thp_limit) {
2048 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
2049 if (thp_update_iterations) {
2050 computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
2051 rates, deferred_logger);
2053 computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
2054 rates, deferred_logger);
2056 auto& ws = well_state.
well(this->name());
2057 ws.surface_rates = rates;
2058 ws.bhp = *bhp_at_thp_limit;
2059 ws.thp = this->getTHPConstraint(summary_state);
2066 template <
typename TypeTag>
2071 const std::vector<Scalar>& mobility,
2074 const int np = this->number_of_phases_;
2075 for (
int p = 0; p < np; ++p) {
2078 const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
2079 const auto connMob =
2080 mobility[FluidSystem::activePhaseToActiveCompIdx(p)] * fs.invB(canonical_phase_idx).value();
2082 connPI[p] = connPICalc(connMob);
2085 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2086 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2088 const auto io = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2089 const auto ig = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2091 const auto vapoil = connPI[ig] * fs.Rv().value();
2092 const auto disgas = connPI[io] * fs.Rs().value();
2094 connPI[io] += vapoil;
2095 connPI[ig] += disgas;
2100 template <
typename TypeTag>
2104 const Phase preferred_phase,
2106 const std::vector<Scalar>& mobility,
2111 if (preferred_phase == Phase::GAS) {
2112 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2114 else if (preferred_phase == Phase::OIL) {
2115 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2117 else if (preferred_phase == Phase::WATER) {
2118 phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
2122 fmt::format(
"Unsupported Injector Type ({}) "
2123 "for well {} during connection I.I. calculation",
2124 static_cast<int>(preferred_phase), this->name()),
2128 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2129 const int canonicalPhaseIdx = FluidSystem::activeToCanonicalPhaseIdx(phase_pos);
2130 connII[phase_pos] = connIICalc(mt * fs.invB(canonicalPhaseIdx).value());
2133 template<
typename TypeTag>
2134 template<
class GasLiftSingleWell>
2135 std::unique_ptr<GasLiftSingleWell>
2144 auto& comm = simulator.vanguard().grid().comm();
2145 ecl_well_map.try_emplace(this->name(), &(this->wellEcl()), this->indexOfWell());
2148 simulator.vanguard().schedule(),
2149 simulator.vanguard().summaryState(),
2150 simulator.episodeIndex(),
2151 simulator.model().newtonMethod().numIterations(),
2160 std::set<int> sync_groups;
2161 const auto& summary_state = simulator.vanguard().summaryState();
2162 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:304
int number_of_local_perforations_
Definition: WellInterfaceGeneric.hpp:340
FluidSystem::Scalar wsolvent_
Definition: WellInterfaceGeneric.hpp:380
const std::vector< FluidSystem::Scalar > & wellIndex() const
Definition: WellInterfaceGeneric.hpp:152
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:701
void checkWellOperability(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1039
Scalar woxygen() const
Definition: WellInterface_impl.hpp:165
IndividualOrGroup
Definition: WellInterface.hpp:246
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:1874
void solveWellEquation(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:805
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:275
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:503
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:136
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:668
std::vector< RateVector > connectionRates_
Definition: WellInterface.hpp:369
void computeConnLevelProdInd(const FluidState &fs, const std::function< Scalar(const Scalar)> &connPICalc, const std::vector< Scalar > &mobility, Scalar *connPI) const
Definition: WellInterface_impl.hpp:2069
void addCellRates(RateVector &rates, int cellIdx) const
Definition: WellInterface_impl.hpp:1003
void gliftBeginTimeStepWellTestUpdateALQ(const Simulator &simulator, WellStateType &well_state, const GroupState< Scalar > &group_state, GLiftEclWells &ecl_well_map, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1066
void assembleWellEq(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:867
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const
Definition: WellInterface_impl.hpp:1019
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:1960
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:883
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:1636
bool updateWellStateWithTHPTargetProd(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:2038
void updateConnectionDFactor(const Simulator &simulator, SingleWellStateType &ws) const
Definition: WellInterface_impl.hpp:1854
Eval getPerfCellPressure(const FluidState &fs) const
Definition: WellInterface_impl.hpp:1945
void initializeProducerWellState(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1686
bool wellUnderZeroGroupRateTarget(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger, std::optional< bool > group_control=std::nullopt) const
Definition: WellInterface_impl.hpp:1604
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:1916
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:2103
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:2137
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:383
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:903
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:1586
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:759
void updateWellOperability(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1128
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:566
Scalar wpolymer() const
Definition: WellInterface_impl.hpp:111
bool updateWellOperabilityFromWellEq(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:1165
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:738
Scalar wmicrobes() const
Definition: WellInterface_impl.hpp:153
virtual void scaleSegmentRatesAndPressure(WellStateType &well_state) const
Definition: WellInterface_impl.hpp:1184
virtual void updateWellStateWithTarget(const Simulator &simulator, const GroupState< Scalar > &group_state, WellStateType &well_state, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1192
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:1623
Definition: WellState.hpp:66
constexpr int numPhases() const
The number of phases present.
Definition: WellState.hpp:249
const SingleWellState< Scalar, IndexTraits > & well(std::size_t well_index) const
Definition: WellState.hpp:290
bool isOpen(const std::string &name) const
Definition: WellState.hpp:197
@ NONE
Definition: DeferredLogger.hpp:46
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilbioeffectsmodules.hh:43
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Static data associated with a well perforation.
Definition: PerforationData.hpp:30