23#ifndef OPM_WELLINTERFACE_HEADER_INCLUDED
25#define OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
26#include <opm/simulators/wells/WellInterface.hpp>
29#include <opm/common/Exceptions.hpp>
31#include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
32#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
34#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
36#include <opm/simulators/wells/GroupState.hpp>
37#include <opm/simulators/wells/TargetCalculator.hpp>
38#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
39#include <opm/simulators/wells/WellHelpers.hpp>
41#include <dune/common/version.hh>
48#include <fmt/format.h>
54 template<
typename TypeTag>
59 const ModelParameters& param,
60 const RateConverterType& rate_converter,
61 const int pvtRegionIdx,
62 const int num_components,
64 const int index_of_well,
77 connectionRates_.resize(this->number_of_perforations_);
79 if constexpr (has_solvent || has_zFraction) {
80 if (well.isInjector()) {
81 auto injectorType = this->well_ecl_.injectorType();
82 if (injectorType == InjectorType::GAS) {
83 this->wsolvent_ = this->well_ecl_.getSolventFraction();
90 template<
typename TypeTag>
94 const std::vector<Scalar>& ,
95 const Scalar gravity_arg,
96 const std::vector<Scalar>& B_avg,
97 const bool changed_to_open_this_step)
99 this->phase_usage_ = phase_usage_arg;
100 this->gravity_ = gravity_arg;
102 this->changed_to_open_this_step_ = changed_to_open_this_step;
108 template<
typename TypeTag>
109 typename WellInterface<TypeTag>::Scalar
110 WellInterface<TypeTag>::
113 if constexpr (has_polymer) {
114 return this->wpolymer_();
124 template<
typename TypeTag>
125 typename WellInterface<TypeTag>::Scalar
126 WellInterface<TypeTag>::
129 if constexpr (has_foam) {
130 return this->wfoam_();
138 template<
typename TypeTag>
139 typename WellInterface<TypeTag>::Scalar
140 WellInterface<TypeTag>::
143 if constexpr (has_brine) {
144 return this->wsalt_();
150 template<
typename TypeTag>
151 typename WellInterface<TypeTag>::Scalar
152 WellInterface<TypeTag>::
155 if constexpr (has_micp) {
156 return this->wmicrobes_();
162 template<
typename TypeTag>
163 typename WellInterface<TypeTag>::Scalar
164 WellInterface<TypeTag>::
167 if constexpr (has_micp) {
168 return this->woxygen_();
180 template<
typename TypeTag>
181 typename WellInterface<TypeTag>::Scalar
182 WellInterface<TypeTag>::
185 if constexpr (has_micp) {
186 return this->wurea_();
192 template<
typename TypeTag>
194 WellInterface<TypeTag>::
195 updateWellControl(
const Simulator& simulator,
196 const IndividualOrGroup iog,
197 WellState<Scalar>& well_state,
198 const GroupState<Scalar>& group_state,
199 DeferredLogger& deferred_logger)
201 if (stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
205 const auto& summaryState = simulator.vanguard().summaryState();
206 const auto& schedule = simulator.vanguard().schedule();
207 const auto& well = this->well_ecl_;
208 auto& ws = well_state.well(this->index_of_well_);
210 if (well.isInjector()) {
211 from = WellInjectorCMode2String(ws.injection_cmode);
213 from = WellProducerCMode2String(ws.production_cmode);
215 bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
219 bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
221 std::ostringstream ss;
222 ss <<
" The control mode for well " << this->name()
223 <<
" is oscillating\n"
224 <<
" We don't allow for more than "
225 << this->param_.max_number_of_well_switches_
226 <<
" switches. The control is kept at " << from;
227 deferred_logger.info(ss.str());
229 this->well_control_log_.push_back(from);
233 bool changed =
false;
234 if (iog == IndividualOrGroup::Individual) {
235 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
236 }
else if (iog == IndividualOrGroup::Group) {
237 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
239 assert(iog == IndividualOrGroup::Both);
240 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
242 Parallel::Communication cc = simulator.vanguard().grid().comm();
246 if (well.isInjector()) {
247 to = WellInjectorCMode2String(ws.injection_cmode);
249 to = WellProducerCMode2String(ws.production_cmode);
251 std::ostringstream ss;
252 ss <<
" Switching control mode for well " << this->name()
256 ss <<
" on rank " << cc.rank();
258 deferred_logger.debug(ss.str());
260 this->well_control_log_.push_back(from);
261 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
262 updatePrimaryVariables(simulator, well_state, deferred_logger);
268 template<
typename TypeTag>
270 WellInterface<TypeTag>::
271 updateWellControlAndStatusLocalIteration(
const Simulator& simulator,
272 WellState<Scalar>& well_state,
273 const GroupState<Scalar>& group_state,
274 const Well::InjectionControls& inj_controls,
275 const Well::ProductionControls& prod_controls,
276 const Scalar wqTotal,
277 DeferredLogger& deferred_logger,
278 const bool fixed_control,
279 const bool fixed_status)
281 const auto& summary_state = simulator.vanguard().summaryState();
282 const auto& schedule = simulator.vanguard().schedule();
283 auto& ws = well_state.well(this->index_of_well_);
285 if (this->isInjector()) {
286 from = WellInjectorCMode2String(ws.injection_cmode);
288 from = WellProducerCMode2String(ws.production_cmode);
290 const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
292 if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
296 const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
297 if (!this->wellIsStopped()){
298 if (wqTotal*sgn <= 0.0 && !fixed_status){
302 bool changed =
false;
303 if (!fixed_control) {
308 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
309 prod_controls.hasControl(Well::ProducerCMode::GRUP);
311 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
312 if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
313 changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
317 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
318 ws.production_cmode == Well::ProducerCMode::THP;
319 if (!thp_controlled){
321 updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
323 ws.thp = this->getTHPConstraint(summary_state);
325 updatePrimaryVariables(simulator, well_state, deferred_logger);
330 }
else if (!fixed_status){
332 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
333 Scalar prod_limit = prod_controls.bhp_limit;
334 Scalar inj_limit = inj_controls.bhp_limit;
335 const bool has_thp = this->wellHasTHPConstraints(summary_state);
337 std::vector<Scalar> rates(this->num_components_);
338 if (this->isInjector()){
339 const Scalar bhp_thp = WellBhpThpCalculator(*this).
340 calculateBhpFromThp(well_state, rates,
343 this->getRefDensity(),
345 inj_limit = std::min(bhp_thp,
static_cast<Scalar
>(inj_controls.bhp_limit));
349 const Scalar bhp_min = WellBhpThpCalculator(*this).
350 calculateMinimumBhpFromThp(well_state,
353 this->getRefDensity());
354 prod_limit = std::max(bhp_min,
static_cast<Scalar
>(prod_controls.bhp_limit));
357 const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
360 well_state.well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
362 well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
373 template<
typename TypeTag>
375 WellInterface<TypeTag>::
376 wellTesting(
const Simulator& simulator,
377 const double simulation_time,
378 WellState<Scalar>& well_state,
379 const GroupState<Scalar>& group_state,
380 WellTestState& well_test_state,
381 DeferredLogger& deferred_logger)
383 deferred_logger.info(
" well " + this->name() +
" is being tested");
385 WellState<Scalar> well_state_copy = well_state;
386 auto& ws = well_state_copy.well(this->indexOfWell());
388 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
389 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
390 updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
391 initPrimaryVariablesEvaluation();
393 if (this->isProducer()) {
394 const auto& schedule = simulator.vanguard().schedule();
395 const auto report_step = simulator.episodeIndex();
396 const auto& glo = schedule.glo(report_step);
398 gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
402 WellTestState welltest_state_temp;
404 bool testWell =
true;
409 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
410 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
412 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
413 deferred_logger.debug(msg);
418 updateWellOperability(simulator, well_state_copy, deferred_logger);
419 if ( !this->isOperableAndSolvable() ) {
420 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
421 deferred_logger.debug(msg);
424 std::vector<Scalar> potentials;
426 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
427 }
catch (
const std::exception& e) {
428 const std::string msg = fmt::format(
"well {}: computeWellPotentials() "
429 "failed during testing for re-opening: ",
430 this->name(), e.what());
431 deferred_logger.info(msg);
434 const int np = well_state_copy.numPhases();
435 for (
int p = 0; p < np; ++p) {
436 ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
438 const bool under_zero_target = this->wellUnderZeroGroupRateTarget(simulator, well_state_copy, deferred_logger);
439 this->updateWellTestState(well_state_copy.well(this->indexOfWell()),
445 this->closeCompletions(welltest_state_temp);
451 if ( welltest_state_temp.num_closed_wells() > 0 ||
452 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
458 if (!welltest_state_temp.well_is_closed(this->name())) {
459 well_test_state.open_well(this->name());
461 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
462 deferred_logger.info(msg);
465 for (
auto& completion : this->well_ecl_.getCompletions()) {
466 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
467 well_test_state.open_completion(this->name(), completion.first);
471 well_state = well_state_copy;
478 template<
typename TypeTag>
480 WellInterface<TypeTag>::
481 iterateWellEquations(
const Simulator& simulator,
483 WellState<Scalar>& well_state,
484 const GroupState<Scalar>& group_state,
485 DeferredLogger& deferred_logger)
487 const auto& summary_state = simulator.vanguard().summaryState();
488 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
489 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
490 bool converged =
false;
493 if (!this->param_.local_well_solver_control_switching_){
494 converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
496 if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN)) {
497 converged = solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
499 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
503 }
catch (NumericalProblem& e ) {
504 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
505 deferred_logger.warning(
"INNER_ITERATION_FAILED", msg);
511 template<
typename TypeTag>
513 WellInterface<TypeTag>::
514 solveWellWithTHPConstraint(
const Simulator& simulator,
516 const Well::InjectionControls& inj_controls,
517 const Well::ProductionControls& prod_controls,
518 WellState<Scalar>& well_state,
519 const GroupState<Scalar>& group_state,
520 DeferredLogger& deferred_logger)
522 const auto& summary_state = simulator.vanguard().summaryState();
523 bool is_operable =
true;
524 bool converged =
true;
525 auto& ws = well_state.well(this->index_of_well_);
527 if (this->wellIsStopped()) {
529 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
530 if (!bhp_target.has_value()) {
532 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
533 deferred_logger.debug(msg);
536 solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
540 ws.thp = this->getTHPConstraint(summary_state);
541 const Scalar bhp = std::max(bhp_target.value(),
542 static_cast<Scalar
>(prod_controls.bhp_limit));
543 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
548 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
551 const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
553 if (converged && !stoppedOrZeroRateTarget(simulator, well_state, deferred_logger) && isThp) {
554 auto rates = well_state.well(this->index_of_well_).surface_rates;
555 this->adaptRatesForVFP(rates);
556 this->updateIPRImplicit(simulator, well_state, deferred_logger);
557 bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state);
560 this->operability_status_.use_vfpexplicit =
true;
561 auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
563 const Scalar reltol = 1e-3;
564 const Scalar cur_bhp = ws.bhp;
565 if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
566 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
567 deferred_logger.debug(msg);
568 solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger);
570 ws.thp = this->getTHPConstraint(summary_state);
571 converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
578 this->operability_status_.use_vfpexplicit =
true;
580 auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger);
581 if (!bhp_target.has_value()) {
585 converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger);
589 const Scalar bhp = std::max(bhp_target.value(),
590 static_cast<Scalar
>(prod_controls.bhp_limit));
591 solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
592 ws.thp = this->getTHPConstraint(summary_state);
593 converged = this->iterateWellEqWithSwitching(simulator, dt,
602 is_operable = is_operable && !this->wellIsStopped();
603 this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
604 this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
608 template<
typename TypeTag>
609 std::optional<typename WellInterface<TypeTag>::Scalar>
610 WellInterface<TypeTag>::
611 estimateOperableBhp(
const Simulator& simulator,
613 WellState<Scalar>& well_state,
614 const SummaryState& summary_state,
615 DeferredLogger& deferred_logger)
619 Scalar bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
621 const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
622 if (!converged || this->wellIsStopped()) {
625 this->updateIPRImplicit(simulator, well_state, deferred_logger);
626 auto rates = well_state.well(this->index_of_well_).surface_rates;
627 this->adaptRatesForVFP(rates);
628 return WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
631 template<
typename TypeTag>
633 WellInterface<TypeTag>::
634 solveWellWithBhp(
const Simulator& simulator,
637 WellState<Scalar>& well_state,
638 DeferredLogger& deferred_logger)
641 auto group_state = GroupState<Scalar>();
642 auto inj_controls = Well::InjectionControls(0);
643 auto prod_controls = Well::ProductionControls(0);
644 auto& ws = well_state.well(this->index_of_well_);
645 auto cmode_inj = ws.injection_cmode;
646 auto cmode_prod = ws.production_cmode;
647 if (this->isInjector()) {
648 inj_controls.addControl(Well::InjectorCMode::BHP);
649 inj_controls.bhp_limit = bhp;
650 inj_controls.cmode = Well::InjectorCMode::BHP;
651 ws.injection_cmode = Well::InjectorCMode::BHP;
653 prod_controls.addControl(Well::ProducerCMode::BHP);
654 prod_controls.bhp_limit = bhp;
655 prod_controls.cmode = Well::ProducerCMode::BHP;
656 ws.production_cmode = Well::ProducerCMode::BHP;
661 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true);
662 ws.injection_cmode = cmode_inj;
663 ws.production_cmode = cmode_prod;
667 template<
typename TypeTag>
669 WellInterface<TypeTag>::
670 solveWellWithZeroRate(
const Simulator& simulator,
672 WellState<Scalar>& well_state,
673 DeferredLogger& deferred_logger)
676 const auto well_status_orig = this->wellStatus_;
679 auto group_state = GroupState<Scalar>();
680 auto inj_controls = Well::InjectionControls(0);
681 auto prod_controls = Well::ProductionControls(0);
682 const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger,
true,
true);
683 this->wellStatus_ = well_status_orig;
687 template<
typename TypeTag>
689 WellInterface<TypeTag>::
690 solveWellForTesting(
const Simulator& simulator,
691 WellState<Scalar>& well_state,
692 const GroupState<Scalar>& group_state,
693 DeferredLogger& deferred_logger)
696 const WellState<Scalar> well_state0 = well_state;
697 const double dt = simulator.timeStepSize();
698 const auto& summary_state = simulator.vanguard().summaryState();
699 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
702 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
703 converged = gliftBeginTimeStepWellTestIterateWellEquations(
704 simulator, dt, well_state, group_state, deferred_logger);
707 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
708 converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
711 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
714 const int max_iter = this->param_.max_welleq_iter_;
715 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
716 + std::to_string(max_iter) +
" iterations");
717 well_state = well_state0;
722 template<
typename TypeTag>
724 WellInterface<TypeTag>::
725 solveWellEquation(
const Simulator& simulator,
726 WellState<Scalar>& well_state,
727 const GroupState<Scalar>& group_state,
728 DeferredLogger& deferred_logger)
730 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
734 const WellState<Scalar> well_state0 = well_state;
735 const double dt = simulator.timeStepSize();
736 bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
746 auto& ws = well_state.well(this->indexOfWell());
747 bool thp_control =
false;
748 if (this->well_ecl_.isInjector()) {
749 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
751 ws.injection_cmode = Well::InjectorCMode::BHP;
752 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
755 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
757 ws.production_cmode = Well::ProducerCMode::BHP;
758 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
762 const std::string msg = std::string(
"The newly opened well ") + this->name()
763 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
764 deferred_logger.debug(msg);
765 converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
770 const int max_iter = this->param_.max_welleq_iter_;
771 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
772 + std::to_string(max_iter) +
" iterations");
773 well_state = well_state0;
779 template <
typename TypeTag>
781 WellInterface<TypeTag>::
782 assembleWellEq(
const Simulator& simulator,
784 WellState<Scalar>& well_state,
785 const GroupState<Scalar>& group_state,
786 DeferredLogger& deferred_logger)
788 prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
789 assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
794 template <
typename TypeTag>
796 WellInterface<TypeTag>::
797 assembleWellEqWithoutIteration(
const Simulator& simulator,
799 WellState<Scalar>& well_state,
800 const GroupState<Scalar>& group_state,
801 DeferredLogger& deferred_logger)
803 const auto& summary_state = simulator.vanguard().summaryState();
804 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
805 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
808 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
813 template<
typename TypeTag>
815 WellInterface<TypeTag>::
816 prepareWellBeforeAssembling(
const Simulator& simulator,
818 WellState<Scalar>& well_state,
819 const GroupState<Scalar>& group_state,
820 DeferredLogger& deferred_logger)
822 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
824 if (this->param_.check_well_operability_iter_)
825 checkWellOperability(simulator, well_state, deferred_logger);
828 const int iteration_idx = simulator.model().newtonMethod().numIterations();
829 if (iteration_idx < this->param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
830 const auto& ws = well_state.well(this->indexOfWell());
831 const auto pmode_orig = ws.production_cmode;
832 const auto imode_orig = ws.injection_cmode;
834 const bool nonzero_rate_original =
835 std::any_of(ws.surface_rates.begin(),
836 ws.surface_rates.begin() + well_state.numPhases(),
837 [](Scalar rate) { return rate != Scalar(0.0); });
839 this->operability_status_.solvable =
true;
840 bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
843 const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
844 if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
848 this->operability_status_.resetOperability();
850 deferred_logger.debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
853 if (ws.production_cmode != pmode_orig || ws.injection_cmode != imode_orig) {
855 if (this->isInjector()) {
856 from = WellInjectorCMode2String(imode_orig);
857 to = WellInjectorCMode2String(ws.injection_cmode);
859 from = WellProducerCMode2String(pmode_orig);
860 to = WellProducerCMode2String(ws.production_cmode);
862 deferred_logger.debug(
" " + this->name() +
" switched from " + from +
" to " + to +
" during local solve");
867 if (this->param_.shut_unsolvable_wells_)
868 this->operability_status_.solvable =
false;
871 if (this->operability_status_.has_negative_potentials) {
872 auto well_state_copy = well_state;
873 std::vector<Scalar> potentials;
875 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
876 }
catch (
const std::exception& e) {
877 const std::string msg = fmt::format(
"well {}: computeWellPotentials() failed "
878 "during attempt to recompute potentials for well: ",
879 this->name(), e.what());
880 deferred_logger.info(msg);
881 this->operability_status_.has_negative_potentials =
true;
883 auto& ws = well_state.well(this->indexOfWell());
884 const int np = well_state.numPhases();
885 for (
int p = 0; p < np; ++p) {
886 ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
889 this->changed_to_open_this_step_ =
false;
890 const bool well_operable = this->operability_status_.isOperableAndSolvable();
892 if (!well_operable && old_well_operable) {
893 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
895 changed_to_stopped_this_step_ =
true;
896 }
else if (well_operable && !old_well_operable) {
897 deferred_logger.info(
" well " + this->name() +
" gets REVIVED during iteration ");
899 changed_to_stopped_this_step_ =
false;
900 this->changed_to_open_this_step_ =
true;
904 template<
typename TypeTag>
906 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
908 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
911 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
912 if (this->cells()[perfIdx] == cellIdx) {
913 for (
int i = 0; i < RateVector::dimension; ++i) {
914 rates[i] += connectionRates_[perfIdx][i];
920 template<
typename TypeTag>
921 typename WellInterface<TypeTag>::Scalar
922 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const
924 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
925 if (this->cells()[perfIdx] == cellIdx) {
926 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
927 return connectionRates_[perfIdx][activeCompIdx].value();
931 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
932 +
" does not perforate cell " + std::to_string(cellIdx));
939 template<
typename TypeTag>
941 WellInterface<TypeTag>::
942 checkWellOperability(
const Simulator& simulator,
943 const WellState<Scalar>& well_state,
944 DeferredLogger& deferred_logger)
947 if (!this->param_.check_well_operability_) {
951 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
955 updateWellOperability(simulator, well_state, deferred_logger);
956 if (!this->operability_status_.isOperableAndSolvable()) {
957 this->operability_status_.use_vfpexplicit =
true;
958 deferred_logger.debug(
"EXPLICIT_LOOKUP_VFP",
959 "well not operable, trying with explicit vfp lookup: " + this->name());
960 updateWellOperability(simulator, well_state, deferred_logger);
964 template<
typename TypeTag>
966 WellInterface<TypeTag>::
967 gliftBeginTimeStepWellTestIterateWellEquations(
const Simulator& simulator,
969 WellState<Scalar>& well_state,
970 const GroupState<Scalar>& group_state,
971 DeferredLogger& deferred_logger)
973 const auto& well_name = this->name();
974 assert(this->wellHasTHPConstraints(simulator.vanguard().summaryState()));
975 const auto& schedule = simulator.vanguard().schedule();
976 auto report_step_idx = simulator.episodeIndex();
977 const auto& glo = schedule.glo(report_step_idx);
978 if(glo.active() && glo.has_well(well_name)) {
979 const auto increment = glo.gaslift_increment();
980 auto alq = well_state.getALQ(well_name);
983 well_state.setALQ(well_name, alq);
985 iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger)))
994 return iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
998 template<
typename TypeTag>
1000 WellInterface<TypeTag>::
1001 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& simulator,
1002 WellState<Scalar>& well_state,
1003 DeferredLogger& deferred_logger)
1005 const auto& summary_state = simulator.vanguard().summaryState();
1006 const auto& well_name = this->name();
1007 if (!this->wellHasTHPConstraints(summary_state)) {
1008 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1009 deferred_logger.info(msg);
1012 const auto& schedule = simulator.vanguard().schedule();
1013 const auto report_step_idx = simulator.episodeIndex();
1014 const auto& glo = schedule.glo(report_step_idx);
1015 if (!glo.has_well(well_name)) {
1016 const std::string msg = fmt::format(
1017 "GLIFT WTEST: Well {} : Gas Lift not activated: "
1018 "WLIFTOPT is probably missing. Skipping.", well_name);
1019 deferred_logger.info(msg);
1022 const auto& gl_well = glo.well(well_name);
1023 auto& max_alq_optional = gl_well.max_rate();
1025 if (max_alq_optional) {
1026 max_alq = *max_alq_optional;
1029 const auto& well_ecl = this->wellEcl();
1030 const auto& controls = well_ecl.productionControls(summary_state);
1031 const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
1032 const auto& alq_values = table.getALQAxis();
1033 max_alq = alq_values.back();
1035 well_state.setALQ(well_name, max_alq);
1036 const std::string msg = fmt::format(
1037 "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
1038 well_name, max_alq);
1039 deferred_logger.info(msg);
1042 template<
typename TypeTag>
1044 WellInterface<TypeTag>::
1045 updateWellOperability(
const Simulator& simulator,
1046 const WellState<Scalar>& well_state,
1047 DeferredLogger& deferred_logger)
1049 if (this->param_.local_well_solver_control_switching_) {
1050 const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger);
1054 deferred_logger.debug(
"Operability check using well equations did not converge for well "
1055 + this->name() +
", reverting to classical approach." );
1058 this->operability_status_.resetOperability();
1060 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1061 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1062 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1063 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1067 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
1068 if (check_thp || bhp_controlled) {
1069 updateIPR(simulator, deferred_logger);
1070 checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
1074 checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger);
1078 template<
typename TypeTag>
1080 WellInterface<TypeTag>::
1081 updateWellOperabilityFromWellEq(
const Simulator& simulator,
1082 const WellState<Scalar>& well_state,
1083 DeferredLogger& deferred_logger)
1086 assert(this->param_.local_well_solver_control_switching_);
1087 this->operability_status_.resetOperability();
1088 WellState<Scalar> well_state_copy = well_state;
1089 const auto& group_state = simulator.problem().wellModel().groupState();
1090 const double dt = simulator.timeStepSize();
1092 bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger);
1096 template<
typename TypeTag>
1098 WellInterface<TypeTag>::
1099 updateWellStateWithTarget(
const Simulator& simulator,
1100 const GroupState<Scalar>& group_state,
1101 WellState<Scalar>& well_state,
1102 DeferredLogger& deferred_logger)
const
1106 const auto& well = this->well_ecl_;
1107 const int well_index = this->index_of_well_;
1108 auto& ws = well_state.well(well_index);
1110 const int np = well_state.numPhases();
1111 const auto& summaryState = simulator.vanguard().summaryState();
1112 const auto& schedule = simulator.vanguard().schedule();
1114 if (this->wellIsStopped()) {
1115 for (
int p = 0; p<np; ++p) {
1116 ws.surface_rates[p] = 0;
1122 if (this->isInjector() )
1124 const auto& controls = well.injectionControls(summaryState);
1126 InjectorType injectorType = controls.injector_type;
1128 switch (injectorType) {
1129 case InjectorType::WATER:
1131 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
1134 case InjectorType::OIL:
1136 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
1139 case InjectorType::GAS:
1141 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
1145 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1148 const auto current = ws.injection_cmode;
1151 case Well::InjectorCMode::RATE:
1153 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1154 if(this->rsRvInj() > 0) {
1155 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1156 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
1157 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1158 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
1160 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 );
1166 case Well::InjectorCMode::RESV:
1168 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1169 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
1170 const Scalar coeff = convert_coeff[phasePos];
1171 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1175 case Well::InjectorCMode::THP:
1177 auto rates = ws.surface_rates;
1178 Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
1182 this->getRefDensity(),
1185 ws.thp = this->getTHPConstraint(summaryState);
1190 Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1191 if (total_rate <= 0.0)
1192 ws.surface_rates = ws.well_potentials;
1196 case Well::InjectorCMode::BHP:
1198 ws.bhp = controls.bhp_limit;
1199 Scalar total_rate = 0.0;
1200 for (
int p = 0; p<np; ++p) {
1201 total_rate += ws.surface_rates[p];
1206 if (total_rate <= 0.0)
1207 ws.surface_rates = ws.well_potentials;
1211 case Well::InjectorCMode::GRUP:
1213 assert(well.isAvailableForGroupControl());
1214 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1215 const Scalar efficiencyFactor = well.getEfficiencyFactor();
1216 std::optional<Scalar> target =
1217 this->getGroupInjectionTargetRate(group,
1226 ws.surface_rates[phasePos] = *target;
1229 case Well::InjectorCMode::CMODE_UNDEFINED:
1231 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
1241 ws.surface_rates[phasePos] = std::max(Scalar{1.e-7}, ws.surface_rates[phasePos]);
1244 ws.bhp = controls.bhp_limit;
1250 const auto current = ws.production_cmode;
1251 const auto& controls = well.productionControls(summaryState);
1253 case Well::ProducerCMode::ORAT:
1255 Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
1258 if (current_rate > 0.0) {
1259 for (
int p = 0; p<np; ++p) {
1260 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1263 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1264 double control_fraction = fractions[pu.phase_pos[Oil]];
1265 if (control_fraction != 0.0) {
1266 for (
int p = 0; p<np; ++p) {
1267 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1273 case Well::ProducerCMode::WRAT:
1275 Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
1278 if (current_rate > 0.0) {
1279 for (
int p = 0; p<np; ++p) {
1280 ws.surface_rates[p] *= controls.water_rate/current_rate;
1283 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1284 const Scalar control_fraction = fractions[pu.phase_pos[Water]];
1285 if (control_fraction != 0.0) {
1286 for (
int p = 0; p<np; ++p) {
1287 ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
1293 case Well::ProducerCMode::GRAT:
1295 Scalar current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
1298 if (current_rate > 0.0) {
1299 for (
int p = 0; p<np; ++p) {
1300 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1303 const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
1304 const Scalar control_fraction = fractions[pu.phase_pos[Gas]];
1305 if (control_fraction != 0.0) {
1306 for (
int p = 0; p<np; ++p) {
1307 ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
1315 case Well::ProducerCMode::LRAT:
1317 Scalar current_rate = - ws.surface_rates[ pu.phase_pos[Water] ]
1318 - ws.surface_rates[ pu.phase_pos[Oil] ];
1321 if (current_rate > 0.0) {
1322 for (
int p = 0; p<np; ++p) {
1323 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1326 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1327 const Scalar control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
1328 if (control_fraction != 0.0) {
1329 for (
int p = 0; p<np; ++p) {
1330 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1336 case Well::ProducerCMode::CRAT:
1338 OPM_DEFLOG_THROW(std::runtime_error,
1339 fmt::format(
"CRAT control not supported, well {}", this->name()),
1342 case Well::ProducerCMode::RESV:
1344 std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1345 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1346 Scalar total_res_rate = 0.0;
1347 for (
int p = 0; p<np; ++p) {
1348 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1350 if (controls.prediction_mode) {
1353 if (total_res_rate > 0.0) {
1354 for (
int p = 0; p<np; ++p) {
1355 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1358 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1359 for (
int p = 0; p<np; ++p) {
1360 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1364 std::vector<Scalar> hrates(this->number_of_phases_,0.);
1365 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1366 hrates[pu.phase_pos[Water]] = controls.water_rate;
1368 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1369 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
1371 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1372 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
1374 std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
1375 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1376 Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1379 if (total_res_rate > 0.0) {
1380 for (
int p = 0; p<np; ++p) {
1381 ws.surface_rates[p] *= target/total_res_rate;
1384 const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1385 for (
int p = 0; p<np; ++p) {
1386 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1392 case Well::ProducerCMode::BHP:
1394 ws.bhp = controls.bhp_limit;
1395 Scalar total_rate = 0.0;
1396 for (
int p = 0; p<np; ++p) {
1397 total_rate -= ws.surface_rates[p];
1402 if (total_rate <= 0.0){
1403 for (
int p = 0; p<np; ++p) {
1404 ws.surface_rates[p] = -ws.well_potentials[p];
1409 case Well::ProducerCMode::THP:
1411 const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger);
1413 if (!update_success) {
1417 auto rates = ws.surface_rates;
1418 this->adaptRatesForVFP(rates);
1419 const Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1420 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1422 ws.thp = this->getTHPConstraint(summaryState);
1426 const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1427 if (total_rate <= 0.0) {
1428 for (
int p = 0; p < this->number_of_phases_; ++p) {
1429 ws.surface_rates[p] = -ws.well_potentials[p];
1435 case Well::ProducerCMode::GRUP:
1437 assert(well.isAvailableForGroupControl());
1438 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1439 const Scalar efficiencyFactor = well.getEfficiencyFactor();
1440 Scalar scale = this->getGroupProductionTargetRate(group,
1450 for (
int p = 0; p<np; ++p) {
1451 ws.surface_rates[p] *= scale;
1453 ws.trivial_target =
false;
1455 ws.trivial_target =
true;
1459 case Well::ProducerCMode::CMODE_UNDEFINED:
1460 case Well::ProducerCMode::NONE:
1462 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1468 ws.bhp = controls.bhp_limit;
1473 template<
typename TypeTag>
1475 WellInterface<TypeTag>::
1476 wellUnderZeroRateTarget(
const Simulator& simulator,
1477 const WellState<Scalar>& well_state,
1478 DeferredLogger& deferred_logger)
const
1481 const bool isGroupControlled = this->wellUnderGroupControl(well_state.well(this->index_of_well_));
1482 if (!isGroupControlled) {
1484 const auto& summaryState = simulator.vanguard().summaryState();
1485 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1487 return this->wellUnderZeroGroupRateTarget(simulator, well_state, deferred_logger, isGroupControlled);
1491 template <
typename TypeTag>
1493 WellInterface<TypeTag>::wellUnderZeroGroupRateTarget(
const Simulator& simulator,
1494 const WellState<Scalar>& well_state,
1495 DeferredLogger& deferred_logger,
1496 const std::optional<bool> group_control)
const
1499 const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.well(this->index_of_well_)));
1500 if (isGroupControlled) {
1501 const auto& summaryState = simulator.vanguard().summaryState();
1502 const auto& group_state = simulator.problem().wellModel().groupState();
1503 const auto& schedule = simulator.vanguard().schedule();
1504 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state, deferred_logger);
1509 template<
typename TypeTag>
1511 WellInterface<TypeTag>::
1512 stoppedOrZeroRateTarget(
const Simulator& simulator,
1513 const WellState<Scalar>& well_state,
1514 DeferredLogger& deferred_logger)
const
1518 return this->wellIsStopped()
1519 || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
1522 template<
typename TypeTag>
1523 std::vector<typename WellInterface<TypeTag>::Scalar>
1524 WellInterface<TypeTag>::
1525 initialWellRateFractions(
const Simulator& simulator,
1526 const WellState<Scalar>& well_state)
const
1528 const int np = this->number_of_phases_;
1529 std::vector<Scalar> scaling_factor(np);
1530 const auto& ws = well_state.well(this->index_of_well_);
1532 Scalar total_potentials = 0.0;
1533 for (
int p = 0; p<np; ++p) {
1534 total_potentials += ws.well_potentials[p];
1536 if (total_potentials > 0) {
1537 for (
int p = 0; p<np; ++p) {
1538 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1540 return scaling_factor;
1544 Scalar total_tw = 0;
1545 const int nperf = this->number_of_perforations_;
1546 for (
int perf = 0; perf < nperf; ++perf) {
1547 total_tw += this->well_index_[perf];
1549 for (
int perf = 0; perf < nperf; ++perf) {
1550 const int cell_idx = this->well_cells_[perf];
1551 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1552 const auto& fs = intQuants.fluidState();
1553 const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
1554 Scalar total_mobility = 0.0;
1555 for (
int p = 0; p < np; ++p) {
1556 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1557 total_mobility += fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value();
1559 for (
int p = 0; p < np; ++p) {
1560 int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
1561 scaling_factor[p] += well_tw_fraction * fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value() / total_mobility;
1564 return scaling_factor;
1569 template <
typename TypeTag>
1578 auto& ws = well_state.well(this->index_of_well_);
1579 int nonzero_rate_index = -1;
1580 const Scalar floating_point_error_epsilon = 1e-14;
1581 for (
int p = 0; p < this->number_of_phases_; ++p) {
1582 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1583 if (nonzero_rate_index == -1) {
1584 nonzero_rate_index = p;
1593 std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
1595 if (nonzero_rate_index == -1) {
1598 for (
int p = 0; p < this->number_of_phases_; ++p) {
1599 ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
1605 const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1606 const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
1607 if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
1608 for (
int p = 0; p < this->number_of_phases_; ++p) {
1609 if (p != nonzero_rate_index) {
1610 const int comp_idx = this->flowPhaseToModelCompIdx(p);
1611 Scalar& rate = ws.surface_rates[p];
1612 rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
1618 template <
typename TypeTag>
1619 std::vector<typename WellInterface<TypeTag>::Scalar>
1622 const IntensiveQuantities& intQuants,
1623 const Scalar trans_mult,
1629 auto wi = std::vector<Scalar>
1630 (this->num_components_, this->well_index_[perf] * trans_mult);
1632 if constexpr (! Indices::gasEnabled) {
1636 const auto& wdfac = this->well_ecl_.getWDFAC();
1638 if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1642 const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
1649 const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
1650 const Scalar Kh = connection.Kh();
1651 const Scalar scaling = 3.141592653589 * Kh * connection.wpimult();
1652 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1654 const Scalar connection_pressure = ws.perf_data.pressure[perf];
1655 const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1656 const Scalar drawdown = cell_pressure - connection_pressure;
1657 const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1658 const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1660 const Scalar b = 2*scaling/wi[gas_comp_idx];
1661 const Scalar c = -2*scaling*mob_g*drawdown;
1663 Scalar consistent_Q = -1.0e20;
1665 const Scalar r2n = b*b + 4*a*c;
1667 const Scalar rn = std::sqrt(r2n);
1668 const Scalar xn1 = (b-rn)*0.5/a;
1672 const Scalar xn2 = (b+rn)*0.5/a;
1673 if (xn2 <= 0 && xn2 > consistent_Q) {
1679 const Scalar r2p = b*b - 4*a*c;
1681 const Scalar rp = std::sqrt(r2p);
1682 const Scalar xp1 = (rp-b)*0.5/a;
1683 if (xp1 > 0 && xp1 < consistent_Q) {
1686 const Scalar xp2 = -(rp+b)*0.5/a;
1687 if (xp2 > 0 && xp2 < consistent_Q) {
1691 wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1696 template <
typename TypeTag>
1698 WellInterface<TypeTag>::
1699 updateConnectionDFactor(
const Simulator& simulator,
1700 SingleWellState<Scalar>& ws)
const
1702 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1706 auto& d_factor = ws.perf_data.connection_d_factor;
1708 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1709 const int cell_idx = this->well_cells_[perf];
1710 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1712 d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
1716 template <
typename TypeTag>
1717 typename WellInterface<TypeTag>::Scalar
1718 WellInterface<TypeTag>::
1719 computeConnectionDFactor(
const int perf,
1720 const IntensiveQuantities& intQuants,
1721 const SingleWellState<Scalar>& ws)
const
1723 auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
1724 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
1728 auto gas_visc = [connection_pressure = ws.perf_data.pressure[perf],
1729 temperature = ws.temperature,
1730 regIdx = this->pvtRegionIdx(), &intQuants]()
1732 const auto rv = getValue(intQuants.fluidState().Rv());
1734 const auto& gasPvt = FluidSystem::gasPvt();
1739 const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
1740 (regIdx, temperature, connection_pressure);
1742 if (! (rv < rv_sat)) {
1743 return gasPvt.saturatedViscosity(regIdx, temperature,
1744 connection_pressure);
1747 return gasPvt.viscosity(regIdx, temperature, connection_pressure,
1748 rv, getValue(intQuants.fluidState().Rvw()));
1751 const auto& connection = this->well_ecl_.getConnections()
1752 [ws.perf_data.ecl_index[perf]];
1754 return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
1758 template <
typename TypeTag>
1760 WellInterface<TypeTag>::
1761 updateConnectionTransmissibilityFactor(
const Simulator& simulator,
1762 SingleWellState<Scalar>& ws)
const
1764 auto connCF = [&connIx = std::as_const(ws.perf_data.ecl_index),
1765 &conns = this->well_ecl_.getConnections()]
1768 return conns[connIx[perf]].CF();
1771 auto& tmult = ws.perf_data.connection_compaction_tmult;
1772 auto& ctf = ws.perf_data.connection_transmissibility_factor;
1774 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
1775 const int cell_idx = this->well_cells_[perf];
1777 const auto& intQuants = simulator.model()
1778 .intensiveQuantities(cell_idx, 0);
1780 tmult[perf] = simulator.problem()
1781 .template wellTransMultiplier<double>(intQuants, cell_idx);
1783 ctf[perf] = connCF(perf) * tmult[perf];
1788 template<
typename TypeTag>
1789 typename WellInterface<TypeTag>::Eval
1790 WellInterface<TypeTag>::getPerfCellPressure(
const typename WellInterface<TypeTag>::FluidState& fs)
const
1792 if constexpr (Indices::oilEnabled) {
1793 return fs.pressure(FluidSystem::oilPhaseIdx);
1794 }
else if constexpr (Indices::gasEnabled) {
1795 return fs.pressure(FluidSystem::gasPhaseIdx);
1797 return fs.pressure(FluidSystem::waterPhaseIdx);
1801 template <
typename TypeTag>
1802 template<
class Value,
class Callback>
1804 WellInterface<TypeTag>::
1805 getMobility(
const Simulator& simulator,
1807 std::vector<Value>& mob,
1808 Callback& extendEval,
1809 [[maybe_unused]] DeferredLogger& deferred_logger)
const
1811 auto relpermArray = []()
1813 if constexpr (std::is_same_v<Value, Scalar>) {
1814 return std::array<Scalar,3>{};
1816 return std::array<Eval,3>{};
1819 const int cell_idx = this->well_cells_[perf];
1820 assert (
int(mob.size()) == this->num_components_);
1821 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1822 const auto& materialLawManager = simulator.problem().materialLawManager();
1826 const int satid = this->saturation_table_number_[perf] - 1;
1827 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1828 if (satid == satid_elem) {
1829 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1830 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1834 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1835 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1837 if constexpr (has_solvent) {
1838 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1841 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1842 auto relativePerms = relpermArray();
1843 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1846 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1849 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1850 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1854 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1855 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1859 if constexpr (has_solvent) {
1860 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1864 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1865 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1866 const auto& connections = this->well_ecl_.getConnections();
1867 const auto& connection = connections[perf_ecl_index];
1868 if (connection.filterCakeActive()) {
1869 for (
auto& val : mob) {
1870 val *= this->inj_fc_multiplier_[perf];
1877 template<
typename TypeTag>
1879 WellInterface<TypeTag>::
1880 updateWellStateWithTHPTargetProd(
const Simulator& simulator,
1881 WellState<Scalar>& well_state,
1882 DeferredLogger& deferred_logger)
const
1884 const auto& summary_state = simulator.vanguard().summaryState();
1886 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1887 simulator, summary_state, this->getALQ(well_state), deferred_logger);
1888 if (bhp_at_thp_limit) {
1889 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
1890 if (thp_update_iterations) {
1891 computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
1892 rates, deferred_logger);
1894 computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
1895 rates, deferred_logger);
1897 auto& ws = well_state.well(this->name());
1898 ws.surface_rates = rates;
1899 ws.bhp = *bhp_at_thp_limit;
1900 ws.thp = this->getTHPConstraint(summary_state);
1907 template <
typename TypeTag>
1909 WellInterface<TypeTag>::
1910 computeConnLevelProdInd(
const FluidState& fs,
1911 const std::function<Scalar(
const Scalar)>& connPICalc,
1912 const std::vector<Scalar>& mobility,
1913 Scalar* connPI)
const
1916 const int np = this->number_of_phases_;
1917 for (
int p = 0; p < np; ++p) {
1920 const auto connMob =
1921 mobility[this->flowPhaseToModelCompIdx(p)]
1922 * fs.invB(this->flowPhaseToModelPhaseIdx(p)).value();
1924 connPI[p] = connPICalc(connMob);
1927 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1928 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1930 const auto io = pu.phase_pos[Oil];
1931 const auto ig = pu.phase_pos[Gas];
1933 const auto vapoil = connPI[ig] * fs.Rv().value();
1934 const auto disgas = connPI[io] * fs.Rs().value();
1936 connPI[io] += vapoil;
1937 connPI[ig] += disgas;
1942 template <
typename TypeTag>
1944 WellInterface<TypeTag>::
1945 computeConnLevelInjInd(
const FluidState& fs,
1946 const Phase preferred_phase,
1947 const std::function<Scalar(
const Scalar)>& connIICalc,
1948 const std::vector<Scalar>& mobility,
1950 DeferredLogger& deferred_logger)
const
1956 if (preferred_phase == Phase::GAS) {
1957 phase_pos = pu.phase_pos[Gas];
1959 else if (preferred_phase == Phase::OIL) {
1960 phase_pos = pu.phase_pos[Oil];
1962 else if (preferred_phase == Phase::WATER) {
1963 phase_pos = pu.phase_pos[Water];
1966 OPM_DEFLOG_THROW(NotImplemented,
1967 fmt::format(
"Unsupported Injector Type ({}) "
1968 "for well {} during connection I.I. calculation",
1969 static_cast<int>(preferred_phase), this->name()),
1973 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1974 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToModelPhaseIdx(phase_pos)).value());
Definition DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition WGState.hpp:31
Definition WellTest.hpp:35
Definition WellInterfaceIndices.hpp:34
Definition WellInterface.hpp:71
WellInterface(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters ¶m, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Constructor.
Definition WellInterface_impl.hpp:56
void updateWellStateRates(const Simulator &simulator, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const
Modify the well_state's rates if there is only one nonzero rate.
Definition WellInterface_impl.hpp:1572
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:62
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
Static data associated with a well perforation.
Definition WellState.hpp:54
Definition BlackoilPhases.hpp:46