22#include <opm/common/Exceptions.hpp>
24#include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
25#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
27#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
29#include <opm/simulators/wells/GroupState.hpp>
30#include <opm/simulators/wells/TargetCalculator.hpp>
31#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
32#include <opm/simulators/wells/WellHelpers.hpp>
34#include <dune/common/version.hh>
39#include <fmt/format.h>
45 template<
typename TypeTag>
51 const RateConverterType& rate_converter,
52 const int pvtRegionIdx,
53 const int num_components,
55 const int index_of_well,
56 const std::vector<PerforationData>& perf_data)
68 connectionRates_.resize(this->number_of_perforations_);
70 if constexpr (has_solvent || has_zFraction) {
71 if (well.isInjector()) {
72 auto injectorType = this->well_ecl_.injectorType();
73 if (injectorType == InjectorType::GAS) {
74 this->wsolvent_ = this->well_ecl_.getSolventFraction();
81 template<
typename TypeTag>
85 const std::vector<double>& ,
86 const double gravity_arg,
88 const std::vector< Scalar >& B_avg,
89 const bool changed_to_open_this_step)
91 this->phase_usage_ = phase_usage_arg;
92 this->gravity_ = gravity_arg;
94 this->changed_to_open_this_step_ = changed_to_open_this_step;
100 template<
typename TypeTag>
102 WellInterface<TypeTag>::
105 if constexpr (has_polymer) {
106 return this->wpolymer_();
116 template<
typename TypeTag>
118 WellInterface<TypeTag>::
121 if constexpr (has_foam) {
122 return this->wfoam_();
130 template<
typename TypeTag>
132 WellInterface<TypeTag>::
135 if constexpr (has_brine) {
136 return this->wsalt_();
142 template<
typename TypeTag>
144 WellInterface<TypeTag>::
147 if constexpr (has_micp) {
148 return this->wmicrobes_();
154 template<
typename TypeTag>
156 WellInterface<TypeTag>::
159 if constexpr (has_micp) {
160 return this->woxygen_();
172 template<
typename TypeTag>
174 WellInterface<TypeTag>::
177 if constexpr (has_micp) {
178 return this->wurea_();
184 template<
typename TypeTag>
186 WellInterface<TypeTag>::
187 updateWellControl(
const Simulator& simulator,
188 const IndividualOrGroup iog,
189 WellState& well_state,
190 const GroupState& group_state,
191 DeferredLogger& deferred_logger)
193 const auto& summary_state = simulator.vanguard().summaryState();
194 if (this->stopppedOrZeroRateTarget(summary_state, well_state)) {
198 const auto& summaryState = simulator.vanguard().summaryState();
199 const auto& schedule = simulator.vanguard().schedule();
200 const auto& well = this->well_ecl_;
201 auto& ws = well_state.well(this->index_of_well_);
203 if (well.isInjector()) {
204 from = WellInjectorCMode2String(ws.injection_cmode);
206 from = WellProducerCMode2String(ws.production_cmode);
208 bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
212 bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
214 std::ostringstream ss;
215 ss <<
" The control mode for well " << this->name()
216 <<
" is oscillating\n"
217 <<
" We don't allow for more than "
218 << param_.max_number_of_well_switches_
219 <<
" switches. The control is kept at " << from;
220 deferred_logger.info(ss.str());
222 this->well_control_log_.push_back(from);
226 bool changed =
false;
227 if (iog == IndividualOrGroup::Individual) {
228 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
229 }
else if (iog == IndividualOrGroup::Group) {
230 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
232 assert(iog == IndividualOrGroup::Both);
233 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
235 Parallel::Communication cc = simulator.vanguard().grid().comm();
239 if (well.isInjector()) {
240 to = WellInjectorCMode2String(ws.injection_cmode);
242 to = WellProducerCMode2String(ws.production_cmode);
244 std::ostringstream ss;
245 ss <<
" Switching control mode for well " << this->name()
249 ss <<
" on rank " << cc.rank();
251 deferred_logger.debug(ss.str());
253 this->well_control_log_.push_back(from);
254 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
255 updatePrimaryVariables(summaryState, well_state, deferred_logger);
261 template<
typename TypeTag>
263 WellInterface<TypeTag>::
264 updateWellControlAndStatusLocalIteration(
const Simulator& simulator,
265 WellState& well_state,
266 const GroupState& group_state,
267 const Well::InjectionControls& inj_controls,
268 const Well::ProductionControls& prod_controls,
269 const double wqTotal,
270 DeferredLogger& deferred_logger,
271 const bool fixed_control,
272 const bool fixed_status)
274 const auto& summary_state = simulator.vanguard().summaryState();
275 const auto& schedule = simulator.vanguard().schedule();
277 if (this->wellUnderZeroRateTarget(summary_state, well_state) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
281 const double sgn = this->isInjector() ? 1.0 : -1.0;
282 if (!this->wellIsStopped()){
283 if (wqTotal*sgn <= 0.0 && !fixed_status){
287 bool changed =
false;
288 if (!fixed_control) {
289 auto& ws = well_state.well(this->index_of_well_);
290 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
291 prod_controls.hasControl(Well::ProducerCMode::GRUP);
293 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
294 if (hasGroupControl) {
295 changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
299 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
300 ws.production_cmode == Well::ProducerCMode::THP;
301 if (!thp_controlled){
303 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
305 ws.thp = this->getTHPConstraint(summary_state);
307 updatePrimaryVariables(summary_state, well_state, deferred_logger);
312 }
else if (!fixed_status){
314 const double bhp = well_state.well(this->index_of_well_).bhp;
315 double prod_limit = prod_controls.bhp_limit;
316 double inj_limit = inj_controls.bhp_limit;
317 const bool has_thp = this->wellHasTHPConstraints(summary_state);
319 std::vector<double> rates(this->num_components_);
320 if (this->isInjector()){
321 const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
322 inj_limit = std::min(bhp_thp, inj_controls.bhp_limit);
325 const double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
326 prod_limit = std::max(bhp_min, prod_controls.bhp_limit);
329 const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
332 well_state.well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
334 well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
345 template<
typename TypeTag>
347 WellInterface<TypeTag>::
348 wellTesting(
const Simulator& simulator,
349 const double simulation_time,
350 WellState& well_state,
351 const GroupState& group_state,
352 WellTestState& well_test_state,
353 DeferredLogger& deferred_logger)
355 deferred_logger.info(
" well " + this->name() +
" is being tested");
357 WellState well_state_copy = well_state;
358 auto& ws = well_state_copy.well(this->indexOfWell());
360 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
361 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
362 const auto& summary_state = simulator.vanguard().summaryState();
363 updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
364 initPrimaryVariablesEvaluation();
366 if (this->isProducer()) {
367 const auto& schedule = simulator.vanguard().schedule();
368 const auto report_step = simulator.episodeIndex();
369 const auto& glo = schedule.glo(report_step);
371 gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
375 WellTestState welltest_state_temp;
377 bool testWell =
true;
382 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
383 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
385 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
386 deferred_logger.debug(msg);
391 updateWellOperability(simulator, well_state_copy, deferred_logger);
392 if ( !this->isOperableAndSolvable() ) {
393 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
394 deferred_logger.debug(msg);
397 std::vector<double> potentials;
399 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
400 }
catch (
const std::exception& e) {
401 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during testing for re-opening: ") + e.what();
402 deferred_logger.info(msg);
405 const int np = well_state_copy.numPhases();
406 for (
int p = 0; p < np; ++p) {
407 ws.well_potentials[p] = std::max(0.0, potentials[p]);
409 this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time,
false, welltest_state_temp, deferred_logger);
410 this->closeCompletions(welltest_state_temp);
416 if ( welltest_state_temp.num_closed_wells() > 0 ||
417 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
423 if (!welltest_state_temp.well_is_closed(this->name())) {
424 well_test_state.open_well(this->name());
426 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
427 deferred_logger.info(msg);
430 for (
auto& completion : this->well_ecl_.getCompletions()) {
431 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
432 well_test_state.open_completion(this->name(), completion.first);
436 well_state = well_state_copy;
443 template<
typename TypeTag>
445 WellInterface<TypeTag>::
446 iterateWellEquations(
const Simulator& simulator,
448 WellState& well_state,
449 const GroupState& group_state,
450 DeferredLogger& deferred_logger)
452 const auto& summary_state = simulator.vanguard().summaryState();
453 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
454 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
455 bool converged =
false;
458 if (!this->param_.local_well_solver_control_switching_){
459 converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
461 if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN)) {
462 converged = solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
464 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
468 }
catch (NumericalProblem& e ) {
469 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
470 deferred_logger.warning(
"INNER_ITERATION_FAILED", msg);
476 template<
typename TypeTag>
478 WellInterface<TypeTag>::
479 solveWellWithTHPConstraint(
const Simulator& simulator,
481 const Well::InjectionControls& inj_controls,
482 const Well::ProductionControls& prod_controls,
483 WellState& well_state,
484 const GroupState& group_state,
485 DeferredLogger& deferred_logger)
487 const auto& summary_state = simulator.vanguard().summaryState();
488 bool is_operable =
true;
489 bool converged =
true;
490 auto& ws = well_state.well(this->index_of_well_);
492 if (this->wellIsStopped()) {
494 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
495 if (!bhp_target.has_value()) {
497 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
498 deferred_logger.debug(msg);
501 solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
505 ws.thp = this->getTHPConstraint(summary_state);
506 const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
507 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
512 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
515 const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
517 if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && isThp) {
518 auto rates = well_state.well(this->index_of_well_).surface_rates;
519 this->adaptRatesForVFP(rates);
520 this->updateIPRImplicit(simulator, well_state, deferred_logger);
521 bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state);
524 this->operability_status_.use_vfpexplicit =
true;
525 auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
527 const double reltol = 1e-3;
528 const double cur_bhp = ws.bhp;
529 if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
530 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
531 deferred_logger.debug(msg);
532 solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger);
534 ws.thp = this->getTHPConstraint(summary_state);
535 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
542 this->operability_status_.use_vfpexplicit =
true;
544 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
545 if (!bhp_target.has_value()) {
549 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
553 const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
554 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
555 ws.thp = this->getTHPConstraint(summary_state);
556 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
560 is_operable = is_operable && !this->wellIsStopped();
561 this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
562 this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
566 template<
typename TypeTag>
567 std::optional<double>
568 WellInterface<TypeTag>::
569 estimateOperableBhp(
const Simulator& simulator,
571 WellState& well_state,
572 const SummaryState& summary_state,
573 DeferredLogger& deferred_logger)
577 double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
579 const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
580 if (!converged || this->wellIsStopped()) {
583 this->updateIPRImplicit(simulator, well_state, deferred_logger);
584 auto rates = well_state.well(this->index_of_well_).surface_rates;
585 this->adaptRatesForVFP(rates);
586 return WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
589 template<
typename TypeTag>
591 WellInterface<TypeTag>::
592 solveWellWithBhp(
const Simulator& simulator,
595 WellState& well_state,
596 DeferredLogger& deferred_logger)
599 auto group_state = GroupState();
600 auto inj_controls = Well::InjectionControls(0);
601 auto prod_controls = Well::ProductionControls(0);
602 auto& ws = well_state.well(this->index_of_well_);
603 auto cmode_inj = ws.injection_cmode;
604 auto cmode_prod = ws.production_cmode;
605 if (this->isInjector()) {
606 inj_controls.addControl(Well::InjectorCMode::BHP);
607 inj_controls.bhp_limit = bhp;
608 inj_controls.cmode = Well::InjectorCMode::BHP;
609 ws.injection_cmode = Well::InjectorCMode::BHP;
611 prod_controls.addControl(Well::ProducerCMode::BHP);
612 prod_controls.bhp_limit = bhp;
613 prod_controls.cmode = Well::ProducerCMode::BHP;
614 ws.production_cmode = Well::ProducerCMode::BHP;
619 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true);
620 ws.injection_cmode = cmode_inj;
621 ws.production_cmode = cmode_prod;
625 template<
typename TypeTag>
627 WellInterface<TypeTag>::
628 solveWellWithZeroRate(
const Simulator& simulator,
630 WellState& well_state,
631 DeferredLogger& deferred_logger)
634 const auto well_status_orig = this->wellStatus_;
637 auto group_state = GroupState();
638 auto inj_controls = Well::InjectionControls(0);
639 auto prod_controls = Well::ProductionControls(0);
640 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true,
true);
641 this->wellStatus_ = well_status_orig;
645 template<
typename TypeTag>
647 WellInterface<TypeTag>::
648 solveWellForTesting(
const Simulator& simulator, WellState& well_state,
const GroupState& group_state,
649 DeferredLogger& deferred_logger)
652 const WellState well_state0 = well_state;
653 const double dt = simulator.timeStepSize();
654 const auto& summary_state = simulator.vanguard().summaryState();
655 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
658 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
659 converged = gliftBeginTimeStepWellTestIterateWellEquations(
660 simulator, dt, well_state, group_state, deferred_logger);
663 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
664 converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
667 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
670 const int max_iter = param_.max_welleq_iter_;
671 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
672 + std::to_string(max_iter) +
" iterations");
673 well_state = well_state0;
678 template<
typename TypeTag>
680 WellInterface<TypeTag>::
681 solveWellEquation(
const Simulator& simulator,
682 WellState& well_state,
683 const GroupState& group_state,
684 DeferredLogger& deferred_logger)
686 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
690 const WellState well_state0 = well_state;
691 const double dt = simulator.timeStepSize();
692 bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
702 auto& ws = well_state.well(this->indexOfWell());
703 bool thp_control =
false;
704 if (this->well_ecl_.isInjector()) {
705 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
707 ws.injection_cmode = Well::InjectorCMode::BHP;
708 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
711 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
713 ws.production_cmode = Well::ProducerCMode::BHP;
714 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
718 const std::string msg = std::string(
"The newly opened well ") + this->name()
719 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
720 deferred_logger.debug(msg);
721 converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
726 const int max_iter = param_.max_welleq_iter_;
727 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
728 + std::to_string(max_iter) +
" iterations");
729 well_state = well_state0;
735 template <
typename TypeTag>
737 WellInterface<TypeTag>::
738 assembleWellEq(
const Simulator& simulator,
740 WellState& well_state,
741 const GroupState& group_state,
742 DeferredLogger& deferred_logger)
745 prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
747 assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
752 template <
typename TypeTag>
754 WellInterface<TypeTag>::
755 assembleWellEqWithoutIteration(
const Simulator& simulator,
757 WellState& well_state,
758 const GroupState& group_state,
759 DeferredLogger& deferred_logger)
761 const auto& summary_state = simulator.vanguard().summaryState();
762 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
763 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
766 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
771 template<
typename TypeTag>
773 WellInterface<TypeTag>::
774 prepareWellBeforeAssembling(
const Simulator& simulator,
776 WellState& well_state,
777 const GroupState& group_state,
778 DeferredLogger& deferred_logger)
780 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
782 if (param_.check_well_operability_iter_)
783 checkWellOperability(simulator, well_state, deferred_logger);
786 const int iteration_idx = simulator.model().newtonMethod().numIterations();
787 if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
788 this->operability_status_.solvable =
true;
789 bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
793 if (param_.shut_unsolvable_wells_)
794 this->operability_status_.solvable =
false;
797 if (this->operability_status_.has_negative_potentials) {
798 auto well_state_copy = well_state;
799 std::vector<double> potentials;
801 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
802 }
catch (
const std::exception& e) {
803 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
804 deferred_logger.info(msg);
805 this->operability_status_.has_negative_potentials =
true;
807 auto& ws = well_state.well(this->indexOfWell());
808 const int np = well_state.numPhases();
809 for (
int p = 0; p < np; ++p) {
810 ws.well_potentials[p] = std::max(0.0, potentials[p]);
813 this->changed_to_open_this_step_ =
false;
814 const bool well_operable = this->operability_status_.isOperableAndSolvable();
816 if (!well_operable && old_well_operable) {
817 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
819 changed_to_stopped_this_step_ =
true;
820 }
else if (well_operable && !old_well_operable) {
821 deferred_logger.info(
" well " + this->name() +
" gets REVIVED during iteration ");
823 changed_to_stopped_this_step_ =
false;
824 this->changed_to_open_this_step_ =
true;
828 template<
typename TypeTag>
830 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
832 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
835 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
836 if (this->cells()[perfIdx] == cellIdx) {
837 for (
int i = 0; i < RateVector::dimension; ++i) {
838 rates[i] += connectionRates_[perfIdx][i];
844 template<
typename TypeTag>
845 typename WellInterface<TypeTag>::Scalar
846 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const {
847 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
848 if (this->cells()[perfIdx] == cellIdx) {
849 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
850 return connectionRates_[perfIdx][activeCompIdx].value();
854 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
855 +
" does not perforate cell " + std::to_string(cellIdx));
862 template<
typename TypeTag>
864 WellInterface<TypeTag>::
865 checkWellOperability(
const Simulator& simulator,
866 const WellState& well_state,
867 DeferredLogger& deferred_logger)
870 if (!param_.check_well_operability_) {
874 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
878 updateWellOperability(simulator, well_state, deferred_logger);
879 if (!this->operability_status_.isOperableAndSolvable()) {
880 this->operability_status_.use_vfpexplicit =
true;
881 deferred_logger.debug(
"EXPLICIT_LOOKUP_VFP",
882 "well not operable, trying with explicit vfp lookup: " + this->name());
883 updateWellOperability(simulator, well_state, deferred_logger);
887 template<
typename TypeTag>
889 WellInterface<TypeTag>::
890 gliftBeginTimeStepWellTestIterateWellEquations(
891 const Simulator& simulator,
893 WellState& well_state,
894 const GroupState &group_state,
895 DeferredLogger& deferred_logger)
897 const auto& well_name = this->name();
898 assert(this->wellHasTHPConstraints(simulator.vanguard().summaryState()));
899 const auto& schedule = simulator.vanguard().schedule();
900 auto report_step_idx = simulator.episodeIndex();
901 const auto& glo = schedule.glo(report_step_idx);
902 if(glo.active() && glo.has_well(well_name)) {
903 const auto increment = glo.gaslift_increment();
904 auto alq = well_state.getALQ(well_name);
907 well_state.setALQ(well_name, alq);
909 iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger)))
918 return iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
922 template<
typename TypeTag>
924 WellInterface<TypeTag>::
925 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& simulator,
926 WellState& well_state,
927 DeferredLogger& deferred_logger)
929 const auto& summary_state = simulator.vanguard().summaryState();
930 const auto& well_name = this->name();
931 if (!this->wellHasTHPConstraints(summary_state)) {
932 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
933 deferred_logger.info(msg);
936 const auto& schedule = simulator.vanguard().schedule();
937 const auto report_step_idx = simulator.episodeIndex();
938 const auto& glo = schedule.glo(report_step_idx);
939 if (!glo.has_well(well_name)) {
940 const std::string msg = fmt::format(
941 "GLIFT WTEST: Well {} : Gas Lift not activated: "
942 "WLIFTOPT is probably missing. Skipping.", well_name);
943 deferred_logger.info(msg);
946 const auto& gl_well = glo.well(well_name);
947 auto& max_alq_optional = gl_well.max_rate();
949 if (max_alq_optional) {
950 max_alq = *max_alq_optional;
953 const auto& well_ecl = this->wellEcl();
954 const auto& controls = well_ecl.productionControls(summary_state);
955 const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
956 const auto& alq_values = table.getALQAxis();
957 max_alq = alq_values.back();
959 well_state.setALQ(well_name, max_alq);
960 const std::string msg = fmt::format(
961 "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
963 deferred_logger.info(msg);
966 template<
typename TypeTag>
968 WellInterface<TypeTag>::
969 updateWellOperability(
const Simulator& simulator,
970 const WellState& well_state,
971 DeferredLogger& deferred_logger)
973 if (this->param_.local_well_solver_control_switching_) {
974 const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger);
978 deferred_logger.debug(
"Operability check using well equations did not converge for well "
979 + this->name() +
", reverting to classical approach." );
982 this->operability_status_.resetOperability();
984 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
985 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
986 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
987 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
991 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
992 if (check_thp || bhp_controlled) {
993 updateIPR(simulator, deferred_logger);
994 checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
998 checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger);
1002 template<
typename TypeTag>
1004 WellInterface<TypeTag>::
1005 updateWellOperabilityFromWellEq(
const Simulator& simulator,
1006 const WellState& well_state,
1007 DeferredLogger& deferred_logger)
1010 assert(this->param_.local_well_solver_control_switching_);
1011 this->operability_status_.resetOperability();
1012 WellState well_state_copy = well_state;
1013 const auto& group_state = simulator.problem().wellModel().groupState();
1014 const double dt = simulator.timeStepSize();
1016 bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger);
1020 template<
typename TypeTag>
1022 WellInterface<TypeTag>::
1023 updateWellStateWithTarget(
const Simulator& simulator,
1024 const GroupState& group_state,
1025 WellState& well_state,
1026 DeferredLogger& deferred_logger)
const
1030 const auto& well = this->well_ecl_;
1031 const int well_index = this->index_of_well_;
1032 auto& ws = well_state.well(well_index);
1034 const int np = well_state.numPhases();
1035 const auto& summaryState = simulator.vanguard().summaryState();
1036 const auto& schedule = simulator.vanguard().schedule();
1038 if (this->wellIsStopped()) {
1039 for (
int p = 0; p<np; ++p) {
1040 ws.surface_rates[p] = 0;
1046 if (this->isInjector() )
1048 const auto& controls = well.injectionControls(summaryState);
1050 InjectorType injectorType = controls.injector_type;
1052 switch (injectorType) {
1053 case InjectorType::WATER:
1055 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
1058 case InjectorType::OIL:
1060 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
1063 case InjectorType::GAS:
1065 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
1069 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1072 const auto current = ws.injection_cmode;
1075 case Well::InjectorCMode::RATE:
1077 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1078 if(this->rsRvInj() > 0) {
1079 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1080 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
1081 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1082 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
1084 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 );
1090 case Well::InjectorCMode::RESV:
1092 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
1093 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
1094 const double coeff = convert_coeff[phasePos];
1095 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1099 case Well::InjectorCMode::THP:
1101 auto rates = ws.surface_rates;
1102 double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
1106 this->getRefDensity(),
1109 ws.thp = this->getTHPConstraint(summaryState);
1114 double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1115 if (total_rate <= 0.0)
1116 ws.surface_rates = ws.well_potentials;
1120 case Well::InjectorCMode::BHP:
1122 ws.bhp = controls.bhp_limit;
1123 double total_rate = 0.0;
1124 for (
int p = 0; p<np; ++p) {
1125 total_rate += ws.surface_rates[p];
1130 if (total_rate <= 0.0)
1131 ws.surface_rates = ws.well_potentials;
1135 case Well::InjectorCMode::GRUP:
1137 assert(well.isAvailableForGroupControl());
1138 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1139 const double efficiencyFactor = well.getEfficiencyFactor();
1140 std::optional<double> target =
1141 this->getGroupInjectionTargetRate(group,
1150 ws.surface_rates[phasePos] = *target;
1153 case Well::InjectorCMode::CMODE_UNDEFINED:
1155 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
1165 ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
1168 ws.bhp = controls.bhp_limit;
1174 const auto current = ws.production_cmode;
1175 const auto& controls = well.productionControls(summaryState);
1177 case Well::ProducerCMode::ORAT:
1179 double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
1182 if (current_rate > 0.0) {
1183 for (
int p = 0; p<np; ++p) {
1184 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1187 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1188 double control_fraction = fractions[pu.phase_pos[Oil]];
1189 if (control_fraction != 0.0) {
1190 for (
int p = 0; p<np; ++p) {
1191 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1197 case Well::ProducerCMode::WRAT:
1199 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
1202 if (current_rate > 0.0) {
1203 for (
int p = 0; p<np; ++p) {
1204 ws.surface_rates[p] *= controls.water_rate/current_rate;
1207 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1208 double control_fraction = fractions[pu.phase_pos[Water]];
1209 if (control_fraction != 0.0) {
1210 for (
int p = 0; p<np; ++p) {
1211 ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
1217 case Well::ProducerCMode::GRAT:
1219 double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
1222 if (current_rate > 0.0) {
1223 for (
int p = 0; p<np; ++p) {
1224 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1227 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1228 double control_fraction = fractions[pu.phase_pos[Gas]];
1229 if (control_fraction != 0.0) {
1230 for (
int p = 0; p<np; ++p) {
1231 ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
1239 case Well::ProducerCMode::LRAT:
1241 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
1242 - ws.surface_rates[ pu.phase_pos[Oil] ];
1245 if (current_rate > 0.0) {
1246 for (
int p = 0; p<np; ++p) {
1247 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1250 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1251 double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
1252 if (control_fraction != 0.0) {
1253 for (
int p = 0; p<np; ++p) {
1254 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1260 case Well::ProducerCMode::CRAT:
1262 OPM_DEFLOG_THROW(std::runtime_error,
1263 fmt::format(
"CRAT control not supported, well {}", this->name()),
1266 case Well::ProducerCMode::RESV:
1268 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
1269 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1270 double total_res_rate = 0.0;
1271 for (
int p = 0; p<np; ++p) {
1272 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1274 if (controls.prediction_mode) {
1277 if (total_res_rate > 0.0) {
1278 for (
int p = 0; p<np; ++p) {
1279 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1282 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1283 for (
int p = 0; p<np; ++p) {
1284 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1288 std::vector<double> hrates(this->number_of_phases_,0.);
1289 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1290 hrates[pu.phase_pos[Water]] = controls.water_rate;
1292 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1293 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
1295 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1296 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
1298 std::vector<double> hrates_resv(this->number_of_phases_,0.);
1299 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1300 double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1303 if (total_res_rate > 0.0) {
1304 for (
int p = 0; p<np; ++p) {
1305 ws.surface_rates[p] *= target/total_res_rate;
1308 const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
1309 for (
int p = 0; p<np; ++p) {
1310 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1317 case Well::ProducerCMode::BHP:
1319 ws.bhp = controls.bhp_limit;
1320 double total_rate = 0.0;
1321 for (
int p = 0; p<np; ++p) {
1322 total_rate -= ws.surface_rates[p];
1327 if (total_rate <= 0.0){
1328 for (
int p = 0; p<np; ++p) {
1329 ws.surface_rates[p] = -ws.well_potentials[p];
1334 case Well::ProducerCMode::THP:
1336 const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger);
1338 if (!update_success) {
1342 auto rates = ws.surface_rates;
1343 this->adaptRatesForVFP(rates);
1344 const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1345 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1347 ws.thp = this->getTHPConstraint(summaryState);
1351 const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1352 if (total_rate <= 0.0) {
1353 for (
int p = 0; p < this->number_of_phases_; ++p) {
1354 ws.surface_rates[p] = -ws.well_potentials[p];
1360 case Well::ProducerCMode::GRUP:
1362 assert(well.isAvailableForGroupControl());
1363 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1364 const double efficiencyFactor = well.getEfficiencyFactor();
1365 double scale = this->getGroupProductionTargetRate(group,
1375 for (
int p = 0; p<np; ++p) {
1376 ws.surface_rates[p] *= scale;
1378 ws.trivial_target =
false;
1380 ws.trivial_target =
true;
1384 case Well::ProducerCMode::CMODE_UNDEFINED:
1385 case Well::ProducerCMode::NONE:
1387 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1394 ws.bhp = controls.bhp_limit;
1399 template<
typename TypeTag>
1401 WellInterface<TypeTag>::
1402 initialWellRateFractions(
const Simulator& simulator,
const WellState& well_state)
const
1404 const int np = this->number_of_phases_;
1405 std::vector<double> scaling_factor(np);
1406 const auto& ws = well_state.well(this->index_of_well_);
1408 double total_potentials = 0.0;
1409 for (
int p = 0; p<np; ++p) {
1410 total_potentials += ws.well_potentials[p];
1412 if (total_potentials > 0) {
1413 for (
int p = 0; p<np; ++p) {
1414 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1416 return scaling_factor;
1420 double total_tw = 0;
1421 const int nperf = this->number_of_perforations_;
1422 for (
int perf = 0; perf < nperf; ++perf) {
1423 total_tw += this->well_index_[perf];
1425 for (
int perf = 0; perf < nperf; ++perf) {
1426 const int cell_idx = this->well_cells_[perf];
1427 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1428 const auto& fs = intQuants.fluidState();
1429 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1430 double total_mobility = 0.0;
1431 for (
int p = 0; p < np; ++p) {
1432 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1433 total_mobility += fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value();
1435 for (
int p = 0; p < np; ++p) {
1436 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1437 scaling_factor[p] += well_tw_fraction * fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value() / total_mobility;
1440 return scaling_factor;
1445 template <
typename TypeTag>
1454 auto& ws = well_state.well(this->index_of_well_);
1455 int nonzero_rate_index = -1;
1456 const double floating_point_error_epsilon = 1e-14;
1457 for (
int p = 0; p < this->number_of_phases_; ++p) {
1458 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1459 if (nonzero_rate_index == -1) {
1460 nonzero_rate_index = p;
1469 std::vector<double> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
1471 if (nonzero_rate_index == -1) {
1474 for (
int p = 0; p < this->number_of_phases_; ++p) {
1475 ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
1481 const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1482 const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
1483 if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
1484 for (
int p = 0; p < this->number_of_phases_; ++p) {
1485 if (p != nonzero_rate_index) {
1486 const int comp_idx = this->flowPhaseToModelCompIdx(p);
1487 double& rate = ws.surface_rates[p];
1488 rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
1494 template <
typename TypeTag>
1498 const IntensiveQuantities& intQuants,
1499 const double trans_mult,
1505 auto wi = std::vector<Scalar>
1506 (this->num_components_, this->well_index_[perf] * trans_mult);
1508 if constexpr (! Indices::gasEnabled) {
1512 const auto& wdfac = this->well_ecl_.getWDFAC();
1514 if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1518 const double d = this->computeConnectionDFactor(perf, intQuants, ws);
1525 const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
1526 const double Kh = connection.Kh();
1527 const double scaling = 3.141592653589 * Kh * connection.wpimult();
1528 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1530 const double connection_pressure = ws.perf_data.pressure[perf];
1531 const double cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1532 const double drawdown = cell_pressure - connection_pressure;
1533 const double invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1534 const double mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1536 const double b = 2*scaling/wi[gas_comp_idx];
1537 const double c = -2*scaling*mob_g*drawdown;
1539 double consistent_Q = -1.0e20;
1541 const double r2n = b*b + 4*a*c;
1543 const double rn = std::sqrt(r2n);
1544 const double xn1 = (b-rn)*0.5/a;
1548 const double xn2 = (b+rn)*0.5/a;
1549 if (xn2 <= 0 && xn2 > consistent_Q) {
1555 const double r2p = b*b - 4*a*c;
1557 const double rp = std::sqrt(r2p);
1558 const double xp1 = (rp-b)*0.5/a;
1559 if (xp1 > 0 && xp1 < consistent_Q) {
1562 const double xp2 = -(rp+b)*0.5/a;
1563 if (xp2 > 0 && xp2 < consistent_Q) {
1567 wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1572 template <
typename TypeTag>
1574 WellInterface<TypeTag>::
1575 updateConnectionDFactor(
const Simulator& simulator, SingleWellState& ws)
const
1577 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1581 auto& d_factor = ws.perf_data.connection_d_factor;
1583 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1584 const int cell_idx = this->well_cells_[perf];
1585 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1587 d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
1591 template <
typename TypeTag>
1593 WellInterface<TypeTag>::
1594 computeConnectionDFactor(
const int perf,
1595 const IntensiveQuantities& intQuants,
1596 const SingleWellState& ws)
const
1598 auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
1599 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
1603 auto gas_visc = [connection_pressure = ws.perf_data.pressure[perf],
1604 temperature = ws.temperature,
1605 regIdx = this->pvtRegionIdx(), &intQuants]()
1607 const auto rv = getValue(intQuants.fluidState().Rv());
1609 const auto& gasPvt = FluidSystem::gasPvt();
1614 const double rv_sat = gasPvt.saturatedOilVaporizationFactor
1615 (regIdx, temperature, connection_pressure);
1617 if (! (rv < rv_sat)) {
1618 return gasPvt.saturatedViscosity(regIdx, temperature,
1619 connection_pressure);
1622 return gasPvt.viscosity(regIdx, temperature, connection_pressure,
1623 rv, getValue(intQuants.fluidState().Rvw()));
1626 const auto& connection = this->well_ecl_.getConnections()
1627 [ws.perf_data.ecl_index[perf]];
1629 return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
1633 template <
typename TypeTag>
1635 WellInterface<TypeTag>::
1636 updateConnectionTransmissibilityFactor(
const Simulator& simulator, SingleWellState& ws)
const
1638 auto connCF = [&connIx = std::as_const(ws.perf_data.ecl_index),
1639 &conns = this->well_ecl_.getConnections()]
1642 return conns[connIx[perf]].CF();
1645 auto& tmult = ws.perf_data.connection_compaction_tmult;
1646 auto& ctf = ws.perf_data.connection_transmissibility_factor;
1648 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1649 const int cell_idx = this->well_cells_[perf];
1651 const auto& intQuants = simulator.model()
1652 .intensiveQuantities(cell_idx, 0);
1654 tmult[perf] = simulator.problem()
1655 .template wellTransMultiplier<double>(intQuants, cell_idx);
1657 ctf[perf] = connCF(perf) * tmult[perf];
1662 template<
typename TypeTag>
1663 typename WellInterface<TypeTag>::Eval
1664 WellInterface<TypeTag>::getPerfCellPressure(
const typename WellInterface<TypeTag>::FluidState& fs)
const
1666 if constexpr (Indices::oilEnabled) {
1667 return fs.pressure(FluidSystem::oilPhaseIdx);
1668 }
else if constexpr (Indices::gasEnabled) {
1669 return fs.pressure(FluidSystem::gasPhaseIdx);
1671 return fs.pressure(FluidSystem::waterPhaseIdx);
1675 template <
typename TypeTag>
1676 template<
class Value,
class Callback>
1678 WellInterface<TypeTag>::
1679 getMobility(
const Simulator& simulator,
1681 std::vector<Value>& mob,
1682 Callback& extendEval,
1683 [[maybe_unused]] DeferredLogger& deferred_logger)
const
1685 auto relpermArray = []()
1687 if constexpr (std::is_same_v<Value, Scalar>) {
1688 return std::array<Scalar,3>{};
1690 return std::array<Eval,3>{};
1693 const int cell_idx = this->well_cells_[perf];
1694 assert (
int(mob.size()) == this->num_components_);
1695 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1696 const auto& materialLawManager = simulator.problem().materialLawManager();
1700 const int satid = this->saturation_table_number_[perf] - 1;
1701 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1702 if (satid == satid_elem) {
1703 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1704 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1708 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1709 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1711 if constexpr (has_solvent) {
1712 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1715 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1716 auto relativePerms = relpermArray();
1717 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1720 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1723 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1724 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1728 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1729 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1733 if constexpr (has_solvent) {
1734 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1738 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1739 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1740 const auto& connections = this->well_ecl_.getConnections();
1741 const auto& connection = connections[perf_ecl_index];
1742 if (connection.filterCakeActive()) {
1743 for (
auto& val : mob) {
1744 val *= this->inj_fc_multiplier_[perf];
1751 template<
typename TypeTag>
1753 WellInterface<TypeTag>::
1754 updateWellStateWithTHPTargetProd(
const Simulator& simulator,
1755 WellState& well_state,
1756 DeferredLogger& deferred_logger)
const
1758 const auto& summary_state = simulator.vanguard().summaryState();
1760 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1761 simulator, summary_state, this->getALQ(well_state), deferred_logger);
1762 if (bhp_at_thp_limit) {
1763 std::vector<double> rates(this->number_of_phases_, 0.0);
1764 if (thp_update_iterations) {
1765 computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
1766 rates, deferred_logger);
1768 computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
1769 rates, deferred_logger);
1771 auto& ws = well_state.well(this->name());
1772 ws.surface_rates = rates;
1773 ws.bhp = *bhp_at_thp_limit;
1774 ws.thp = this->getTHPConstraint(summary_state);
1781 template <
typename TypeTag>
1783 WellInterface<TypeTag>::
1784 computeConnLevelProdInd(
const FluidState& fs,
1785 const std::function<
double(
const double)>& connPICalc,
1786 const std::vector<Scalar>& mobility,
1787 double* connPI)
const
1790 const int np = this->number_of_phases_;
1791 for (
int p = 0; p < np; ++p) {
1794 const auto connMob =
1795 mobility[this->flowPhaseToModelCompIdx(p)]
1796 * fs.invB(this->flowPhaseToModelPhaseIdx(p)).value();
1798 connPI[p] = connPICalc(connMob);
1801 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1802 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1804 const auto io = pu.phase_pos[Oil];
1805 const auto ig = pu.phase_pos[Gas];
1807 const auto vapoil = connPI[ig] * fs.Rv().value();
1808 const auto disgas = connPI[io] * fs.Rs().value();
1810 connPI[io] += vapoil;
1811 connPI[ig] += disgas;
1816 template <
typename TypeTag>
1818 WellInterface<TypeTag>::
1819 computeConnLevelInjInd(
const FluidState& fs,
1820 const Phase preferred_phase,
1821 const std::function<
double(
const double)>& connIICalc,
1822 const std::vector<Scalar>& mobility,
1824 DeferredLogger& deferred_logger)
const
1830 if (preferred_phase == Phase::GAS) {
1831 phase_pos = pu.phase_pos[Gas];
1833 else if (preferred_phase == Phase::OIL) {
1834 phase_pos = pu.phase_pos[Oil];
1836 else if (preferred_phase == Phase::WATER) {
1837 phase_pos = pu.phase_pos[Water];
1840 OPM_DEFLOG_THROW(NotImplemented,
1841 fmt::format(
"Unsupported Injector Type ({}) "
1842 "for well {} during connection I.I. calculation",
1843 static_cast<int>(preferred_phase), this->name()),
1847 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1848 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToModelPhaseIdx(phase_pos)).value());
Definition DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition ParallelWellInfo.hpp:184
Definition SingleWellState.hpp:40
Definition WellInterfaceIndices.hpp:35
Definition WellInterface.hpp:75
WellInterface(const Well &well, const ParallelWellInfo &pw_info, const int time_step, const ModelParameters ¶m, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData > &perf_data)
Constructor.
Definition WellInterface_impl.hpp:47
void updateWellStateRates(const Simulator &simulator, WellState &well_state, DeferredLogger &deferred_logger) const
Modify the well_state's rates if there is only one nonzero rate.
Definition WellInterface_impl.hpp:1448
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:61
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
Solver parameters for the BlackoilModel.
Definition BlackoilModelParameters.hpp:484
Definition BlackoilPhases.hpp:46