22#include <opm/common/Exceptions.hpp>
24#include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
25#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
26#include <opm/simulators/wells/GroupState.hpp>
27#include <opm/simulators/wells/TargetCalculator.hpp>
28#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
29#include <opm/simulators/wells/WellHelpers.hpp>
31#include <dune/common/version.hh>
33#include <fmt/format.h>
41 template<
typename TypeTag>
48 const int pvtRegionIdx,
52 const std::vector<PerforationData>& perf_data)
64 connectionRates_.resize(this->number_of_perforations_);
66 if constexpr (has_solvent || has_zFraction) {
67 if (well.isInjector()) {
70 this->wsolvent_ = this->well_ecl_.getSolventFraction();
77 template<
typename TypeTag>
81 const std::vector<double>& ,
84 const std::vector< Scalar >&
B_avg,
96 template<
typename TypeTag>
98 WellInterface<TypeTag>::
101 if constexpr (has_polymer) {
102 return this->wpolymer_();
112 template<
typename TypeTag>
114 WellInterface<TypeTag>::
117 if constexpr (has_foam) {
118 return this->wfoam_();
126 template<
typename TypeTag>
128 WellInterface<TypeTag>::
131 if constexpr (has_brine) {
132 return this->wsalt_();
138 template<
typename TypeTag>
140 WellInterface<TypeTag>::
143 if constexpr (has_micp) {
144 return this->wmicrobes_();
150 template<
typename TypeTag>
152 WellInterface<TypeTag>::
155 if constexpr (has_micp) {
156 return this->woxygen_();
168 template<
typename TypeTag>
170 WellInterface<TypeTag>::
173 if constexpr (has_micp) {
174 return this->wurea_();
180 template<
typename TypeTag>
182 WellInterface<TypeTag>::
183 updateWellControl(
const Simulator& ebos_simulator,
184 const IndividualOrGroup iog,
185 WellState& well_state,
186 const GroupState& group_state,
187 DeferredLogger& deferred_logger)
189 const auto& summary_state = ebos_simulator.vanguard().summaryState();
190 if (this->stopppedOrZeroRateTarget(summary_state, well_state)) {
194 const auto& summaryState = ebos_simulator.vanguard().summaryState();
195 const auto& schedule = ebos_simulator.vanguard().schedule();
196 const auto& well = this->well_ecl_;
197 auto& ws = well_state.well(this->index_of_well_);
199 if (well.isInjector()) {
200 from = WellInjectorCMode2String(ws.injection_cmode);
202 from = WellProducerCMode2String(ws.production_cmode);
204 bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
208 bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
210 std::ostringstream ss;
211 ss <<
" The control model for well " << this->name()
212 <<
" is oscillating\n"
213 <<
" We don't allow for more than "
214 << param_.max_number_of_well_switches_
215 <<
" switches. The control is kept at " << from;
216 deferred_logger.info(ss.str());
218 this->well_control_log_.push_back(from);
222 bool changed =
false;
223 if (iog == IndividualOrGroup::Individual) {
224 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
225 }
else if (iog == IndividualOrGroup::Group) {
226 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
228 assert(iog == IndividualOrGroup::Both);
229 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
231 Parallel::Communication cc = ebos_simulator.vanguard().grid().comm();
235 if (well.isInjector()) {
236 to = WellInjectorCMode2String(ws.injection_cmode);
238 to = WellProducerCMode2String(ws.production_cmode);
240 std::ostringstream ss;
241 ss <<
" Switching control mode for well " << this->name()
245 ss <<
" on rank " << cc.rank();
247 deferred_logger.debug(ss.str());
249 this->well_control_log_.push_back(from);
250 updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
251 updatePrimaryVariables(summaryState, well_state, deferred_logger);
257 template<
typename TypeTag>
259 WellInterface<TypeTag>::
260 updateWellControlAndStatusLocalIteration(
const Simulator& ebos_simulator,
261 WellState& well_state,
262 const GroupState& group_state,
263 const Well::InjectionControls& inj_controls,
264 const Well::ProductionControls& prod_controls,
265 const double wqTotal,
266 DeferredLogger& deferred_logger)
268 const auto& summary_state = ebos_simulator.vanguard().summaryState();
269 const auto& schedule = ebos_simulator.vanguard().schedule();
271 if (this->wellUnderZeroRateTarget(summary_state, well_state) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
275 const double sgn = this->isInjector() ? 1.0 : -1.0;
276 if (!this->wellIsStopped()){
277 if (wqTotal*sgn <= 0.0){
281 bool changed =
false;
282 auto& ws = well_state.well(this->index_of_well_);
283 const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
284 prod_controls.hasControl(Well::ProducerCMode::GRUP);
286 changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
288 if (hasGroupControl) {
289 changed = this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
293 const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
294 ws.production_cmode == Well::ProducerCMode::THP;
295 if (!thp_controlled){
297 updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
299 ws.thp = this->getTHPConstraint(summary_state);
301 updatePrimaryVariables(summary_state, well_state, deferred_logger);
307 const double bhp = well_state.well(this->index_of_well_).bhp;
308 double prod_limit = prod_controls.bhp_limit;
309 double inj_limit = inj_controls.bhp_limit;
310 const bool has_thp = this->wellHasTHPConstraints(summary_state);
316 std::vector<double> rates(this->num_components_);
317 const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
318 if (this->isInjector()){
319 inj_limit = std::min(bhp_thp, inj_controls.bhp_limit);
321 prod_limit = std::max(bhp_thp, prod_controls.bhp_limit);
324 const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
334 template<
typename TypeTag>
336 WellInterface<TypeTag>::
337 wellTesting(
const Simulator& simulator,
338 const double simulation_time,
339 WellState& well_state,
340 const GroupState& group_state,
341 WellTestState& well_test_state,
342 DeferredLogger& deferred_logger)
344 deferred_logger.info(
" well " + this->name() +
" is being tested");
346 WellState well_state_copy = well_state;
347 auto& ws = well_state_copy.well(this->indexOfWell());
349 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
350 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
351 const auto& summary_state = simulator.vanguard().summaryState();
352 updatePrimaryVariables(summary_state, well_state_copy, deferred_logger);
353 initPrimaryVariablesEvaluation();
355 if (this->isProducer()) {
356 gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
359 WellTestState welltest_state_temp;
361 bool testWell =
true;
366 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
367 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
369 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
370 deferred_logger.debug(msg);
375 updateWellOperability(simulator, well_state_copy, deferred_logger);
376 if ( !this->isOperableAndSolvable() ) {
377 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
378 deferred_logger.debug(msg);
381 std::vector<double> potentials;
383 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
384 }
catch (
const std::exception& e) {
385 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during testing for re-opening: ") + e.what();
386 deferred_logger.info(msg);
389 const int np = well_state_copy.numPhases();
390 for (
int p = 0; p < np; ++p) {
391 ws.well_potentials[p] = std::max(0.0, potentials[p]);
393 this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time,
false, welltest_state_temp, deferred_logger);
394 this->closeCompletions(welltest_state_temp);
400 if ( welltest_state_temp.num_closed_wells() > 0 ||
401 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
407 if (!welltest_state_temp.well_is_closed(this->name())) {
408 well_test_state.open_well(this->name());
410 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
411 deferred_logger.info(msg);
414 for (
auto& completion : this->well_ecl_.getCompletions()) {
415 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
416 well_test_state.open_completion(this->name(), completion.first);
420 well_state = well_state_copy;
427 template<
typename TypeTag>
429 WellInterface<TypeTag>::
430 iterateWellEquations(
const Simulator& ebosSimulator,
432 WellState& well_state,
433 const GroupState& group_state,
434 DeferredLogger& deferred_logger)
436 const auto& summary_state = ebosSimulator.vanguard().summaryState();
437 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
438 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
439 bool converged =
false;
442 if (!this->param_.local_well_solver_control_switching_){
443 converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
445 converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
448 }
catch (NumericalProblem& e ) {
449 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
450 deferred_logger.warning(
"INNER_ITERATION_FAILED", msg);
457 template<
typename TypeTag>
459 WellInterface<TypeTag>::
460 solveWellForTesting(
const Simulator& ebosSimulator, WellState& well_state,
const GroupState& group_state,
461 DeferredLogger& deferred_logger)
464 const WellState well_state0 = well_state;
465 const double dt = ebosSimulator.timeStepSize();
466 const auto& summary_state = ebosSimulator.vanguard().summaryState();
467 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
470 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
471 converged = gliftBeginTimeStepWellTestIterateWellEquations(
472 ebosSimulator, dt, well_state, group_state, deferred_logger);
475 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
476 converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
479 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
482 const int max_iter = param_.max_welleq_iter_;
483 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
484 + std::to_string(max_iter) +
" iterations");
485 well_state = well_state0;
490 template<
typename TypeTag>
492 WellInterface<TypeTag>::
493 solveWellEquation(
const Simulator& ebosSimulator,
494 WellState& well_state,
495 const GroupState& group_state,
496 DeferredLogger& deferred_logger)
498 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
502 const WellState well_state0 = well_state;
503 const double dt = ebosSimulator.timeStepSize();
504 bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
514 auto& ws = well_state.well(this->indexOfWell());
515 bool thp_control =
false;
516 if (this->well_ecl_.isInjector()) {
517 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
519 ws.injection_cmode = Well::InjectorCMode::BHP;
520 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
523 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
525 ws.production_cmode = Well::ProducerCMode::BHP;
526 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
530 const std::string msg = std::string(
"The newly opened well ") + this->name()
531 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
532 deferred_logger.debug(msg);
533 converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
538 const int max_iter = param_.max_welleq_iter_;
539 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
540 + std::to_string(max_iter) +
" iterations");
541 well_state = well_state0;
547 template <
typename TypeTag>
549 WellInterface<TypeTag>::
550 assembleWellEq(
const Simulator& ebosSimulator,
552 WellState& well_state,
553 const GroupState& group_state,
554 DeferredLogger& deferred_logger)
557 prepareWellBeforeAssembling(ebosSimulator, dt, well_state, group_state, deferred_logger);
559 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, group_state, deferred_logger);
564 template <
typename TypeTag>
566 WellInterface<TypeTag>::
567 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
569 WellState& well_state,
570 const GroupState& group_state,
571 DeferredLogger& deferred_logger)
573 const auto& summary_state = ebosSimulator.vanguard().summaryState();
574 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
575 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
578 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
583 template<
typename TypeTag>
585 WellInterface<TypeTag>::
586 prepareWellBeforeAssembling(
const Simulator& ebosSimulator,
588 WellState& well_state,
589 const GroupState& group_state,
590 DeferredLogger& deferred_logger)
592 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
594 if (param_.check_well_operability_iter_)
595 checkWellOperability(ebosSimulator, well_state, deferred_logger);
598 const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
599 if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
600 this->operability_status_.solvable =
true;
601 bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
605 if (param_.shut_unsolvable_wells_)
606 this->operability_status_.solvable =
false;
609 if (this->operability_status_.has_negative_potentials) {
610 auto well_state_copy = well_state;
611 std::vector<double> potentials;
613 computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
614 }
catch (
const std::exception& e) {
615 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
616 deferred_logger.info(msg);
617 this->operability_status_.has_negative_potentials =
true;
619 auto& ws = well_state.well(this->indexOfWell());
620 const int np = well_state.numPhases();
621 for (
int p = 0; p < np; ++p) {
622 ws.well_potentials[p] = std::max(0.0, potentials[p]);
625 this->changed_to_open_this_step_ =
false;
626 const bool well_operable = this->operability_status_.isOperableAndSolvable();
628 if (!well_operable && old_well_operable) {
629 if (this->param_.local_well_solver_control_switching_) {
630 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
632 changed_to_stopped_this_step_ =
true;
636 if (this->well_ecl_.getAutomaticShutIn()) {
637 deferred_logger.info(
" well " + this->name() +
" gets SHUT during iteration ");
639 if (!this->wellIsStopped()) {
640 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
642 changed_to_stopped_this_step_ =
true;
646 }
else if (well_operable && !old_well_operable) {
647 deferred_logger.info(
" well " + this->name() +
" gets REVIVED during iteration ");
649 changed_to_stopped_this_step_ =
false;
650 this->changed_to_open_this_step_ =
true;
654 template<
typename TypeTag>
656 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
658 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
661 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
662 if (this->cells()[perfIdx] == cellIdx) {
663 for (
int i = 0; i < RateVector::dimension; ++i) {
664 rates[i] += connectionRates_[perfIdx][i];
670 template<
typename TypeTag>
671 typename WellInterface<TypeTag>::Scalar
672 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const {
673 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
674 if (this->cells()[perfIdx] == cellIdx) {
675 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
676 return connectionRates_[perfIdx][activeCompIdx].value();
680 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
681 +
" does not perforate cell " + std::to_string(cellIdx));
688 template<
typename TypeTag>
690 WellInterface<TypeTag>::
691 checkWellOperability(
const Simulator& ebos_simulator,
692 const WellState& well_state,
693 DeferredLogger& deferred_logger)
696 if (!param_.check_well_operability_) {
700 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
704 updateWellOperability(ebos_simulator, well_state, deferred_logger);
705 if (!this->operability_status_.isOperableAndSolvable()) {
706 this->operability_status_.use_vfpexplicit =
true;
707 deferred_logger.debug(
"EXPLICIT_LOOKUP_VFP",
708 "well not operable, trying with explicit vfp lookup: " + this->name());
709 updateWellOperability(ebos_simulator, well_state, deferred_logger);
713 template<
typename TypeTag>
715 WellInterface<TypeTag>::
716 gliftBeginTimeStepWellTestIterateWellEquations(
717 const Simulator& ebos_simulator,
719 WellState& well_state,
720 const GroupState &group_state,
721 DeferredLogger& deferred_logger)
723 const auto& well_name = this->name();
724 assert(this->wellHasTHPConstraints(ebos_simulator.vanguard().summaryState()));
725 const auto& schedule = ebos_simulator.vanguard().schedule();
726 auto report_step_idx = ebos_simulator.episodeIndex();
727 const auto& glo = schedule.glo(report_step_idx);
728 if(glo.has_well(well_name)) {
729 auto increment = glo.gaslift_increment();
730 auto alq = well_state.getALQ(well_name);
733 well_state.setALQ(well_name, alq);
735 iterateWellEquations(ebos_simulator, dt, well_state, group_state, deferred_logger)))
744 return iterateWellEquations(ebos_simulator, dt, well_state, group_state, deferred_logger);
748 template<
typename TypeTag>
750 WellInterface<TypeTag>::
751 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& ebos_simulator,
752 WellState& well_state,
753 DeferredLogger& deferred_logger)
755 const auto& summary_state = ebos_simulator.vanguard().summaryState();
756 const auto& well_name = this->name();
757 if (!this->wellHasTHPConstraints(summary_state)) {
758 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
759 deferred_logger.info(msg);
762 const auto& well_ecl = this->wellEcl();
763 const auto& schedule = ebos_simulator.vanguard().schedule();
764 auto report_step_idx = ebos_simulator.episodeIndex();
765 const auto& glo = schedule.glo(report_step_idx);
766 if (!glo.has_well(well_name)) {
767 const std::string msg = fmt::format(
768 "GLIFT WTEST: Well {} : Gas Lift not activated: "
769 "WLIFTOPT is probably missing. Skipping.", well_name);
770 deferred_logger.info(msg);
773 const auto& gl_well = glo.well(well_name);
774 auto& max_alq_optional = gl_well.max_rate();
776 if (max_alq_optional) {
777 max_alq = *max_alq_optional;
780 const auto& controls = well_ecl.productionControls(summary_state);
781 const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
782 const auto& alq_values = table.getALQAxis();
783 max_alq = alq_values.back();
785 well_state.setALQ(well_name, max_alq);
786 const std::string msg = fmt::format(
787 "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
789 deferred_logger.info(msg);
792 template<
typename TypeTag>
794 WellInterface<TypeTag>::
795 updateWellOperability(
const Simulator& ebos_simulator,
796 const WellState& well_state,
797 DeferredLogger& deferred_logger)
799 if (this->param_.local_well_solver_control_switching_) {
800 const bool success = updateWellOperabilityFromWellEq(ebos_simulator, well_state, deferred_logger);
804 deferred_logger.debug(
"Operability check using well equations did not converge for well "
805 + this->name() +
", reverting to classical approach." );
808 this->operability_status_.resetOperability();
810 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
811 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
812 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
813 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
817 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
818 if (check_thp || bhp_controlled) {
819 updateIPR(ebos_simulator, deferred_logger);
820 checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
824 checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
828 template<
typename TypeTag>
830 WellInterface<TypeTag>::
831 updateWellOperabilityFromWellEq(
const Simulator& ebos_simulator,
832 const WellState& well_state,
833 DeferredLogger& deferred_logger)
836 assert(this->param_.local_well_solver_control_switching_);
837 this->operability_status_.resetOperability();
838 WellState well_state_copy = well_state;
839 const auto& group_state = ebos_simulator.problem().wellModel().groupState();
840 const double dt = ebos_simulator.timeStepSize();
842 bool converged = iterateWellEquations(ebos_simulator, dt, well_state_copy, group_state, deferred_logger);
846 template<
typename TypeTag>
848 WellInterface<TypeTag>::
849 updateWellStateWithTarget(
const Simulator& ebos_simulator,
850 const GroupState& group_state,
851 WellState& well_state,
852 DeferredLogger& deferred_logger)
const
856 const auto& well = this->well_ecl_;
857 const int well_index = this->index_of_well_;
858 auto& ws = well_state.well(well_index);
860 const int np = well_state.numPhases();
861 const auto& summaryState = ebos_simulator.vanguard().summaryState();
862 const auto& schedule = ebos_simulator.vanguard().schedule();
864 if (this->wellIsStopped()) {
865 for (
int p = 0; p<np; ++p) {
866 ws.surface_rates[p] = 0;
872 if (this->isInjector() )
874 const auto& controls = well.injectionControls(summaryState);
876 InjectorType injectorType = controls.injector_type;
878 switch (injectorType) {
879 case InjectorType::WATER:
881 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
884 case InjectorType::OIL:
886 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
889 case InjectorType::GAS:
891 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
895 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
898 const auto current = ws.injection_cmode;
901 case Well::InjectorCMode::RATE:
903 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
904 if(this->rsRvInj() > 0) {
905 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
906 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
907 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
908 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
910 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 );
916 case Well::InjectorCMode::RESV:
918 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
919 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
920 const double coeff = convert_coeff[phasePos];
921 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
925 case Well::InjectorCMode::THP:
927 auto rates = ws.surface_rates;
928 double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
932 this->getRefDensity(),
935 ws.thp = this->getTHPConstraint(summaryState);
940 double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
941 if (total_rate <= 0.0)
942 ws.surface_rates = ws.well_potentials;
946 case Well::InjectorCMode::BHP:
948 ws.bhp = controls.bhp_limit;
949 double total_rate = 0.0;
950 for (
int p = 0; p<np; ++p) {
951 total_rate += ws.surface_rates[p];
956 if (total_rate <= 0.0)
957 ws.surface_rates = ws.well_potentials;
961 case Well::InjectorCMode::GRUP:
963 assert(well.isAvailableForGroupControl());
964 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
965 const double efficiencyFactor = well.getEfficiencyFactor();
966 std::optional<double> target =
967 this->getGroupInjectionTargetRate(group,
976 ws.surface_rates[phasePos] = *target;
979 case Well::InjectorCMode::CMODE_UNDEFINED:
981 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
991 ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
994 ws.bhp = controls.bhp_limit;
1000 const auto current = ws.production_cmode;
1001 const auto& controls = well.productionControls(summaryState);
1003 case Well::ProducerCMode::ORAT:
1005 double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
1008 if (current_rate > 0.0) {
1009 for (
int p = 0; p<np; ++p) {
1010 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1013 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1014 double control_fraction = fractions[pu.phase_pos[Oil]];
1015 if (control_fraction != 0.0) {
1016 for (
int p = 0; p<np; ++p) {
1017 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1023 case Well::ProducerCMode::WRAT:
1025 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
1028 if (current_rate > 0.0) {
1029 for (
int p = 0; p<np; ++p) {
1030 ws.surface_rates[p] *= controls.water_rate/current_rate;
1033 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1034 double control_fraction = fractions[pu.phase_pos[Water]];
1035 if (control_fraction != 0.0) {
1036 for (
int p = 0; p<np; ++p) {
1037 ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
1043 case Well::ProducerCMode::GRAT:
1045 double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
1048 if (current_rate > 0.0) {
1049 for (
int p = 0; p<np; ++p) {
1050 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1053 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1054 double control_fraction = fractions[pu.phase_pos[Gas]];
1055 if (control_fraction != 0.0) {
1056 for (
int p = 0; p<np; ++p) {
1057 ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
1065 case Well::ProducerCMode::LRAT:
1067 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
1068 - ws.surface_rates[ pu.phase_pos[Oil] ];
1071 if (current_rate > 0.0) {
1072 for (
int p = 0; p<np; ++p) {
1073 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1076 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1077 double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
1078 if (control_fraction != 0.0) {
1079 for (
int p = 0; p<np; ++p) {
1080 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1086 case Well::ProducerCMode::CRAT:
1088 OPM_DEFLOG_THROW(std::runtime_error,
1089 fmt::format(
"CRAT control not supported, well {}", this->name()),
1092 case Well::ProducerCMode::RESV:
1094 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
1095 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1096 double total_res_rate = 0.0;
1097 for (
int p = 0; p<np; ++p) {
1098 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1100 if (controls.prediction_mode) {
1103 if (total_res_rate > 0.0) {
1104 for (
int p = 0; p<np; ++p) {
1105 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1108 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1109 for (
int p = 0; p<np; ++p) {
1110 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1114 std::vector<double> hrates(this->number_of_phases_,0.);
1115 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1116 hrates[pu.phase_pos[Water]] = controls.water_rate;
1118 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1119 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
1121 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1122 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
1124 std::vector<double> hrates_resv(this->number_of_phases_,0.);
1125 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1126 double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1129 if (total_res_rate > 0.0) {
1130 for (
int p = 0; p<np; ++p) {
1131 ws.surface_rates[p] *= target/total_res_rate;
1134 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1135 for (
int p = 0; p<np; ++p) {
1136 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1143 case Well::ProducerCMode::BHP:
1145 ws.bhp = controls.bhp_limit;
1146 double total_rate = 0.0;
1147 for (
int p = 0; p<np; ++p) {
1148 total_rate -= ws.surface_rates[p];
1153 if (total_rate <= 0.0){
1154 for (
int p = 0; p<np; ++p) {
1155 ws.surface_rates[p] = -ws.well_potentials[p];
1160 case Well::ProducerCMode::THP:
1162 const bool update_success = updateWellStateWithTHPTargetProd(ebos_simulator, well_state, deferred_logger);
1164 if (!update_success) {
1168 auto rates = ws.surface_rates;
1169 this->adaptRatesForVFP(rates);
1170 const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1171 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1173 ws.thp = this->getTHPConstraint(summaryState);
1177 const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1178 if (total_rate <= 0.0) {
1179 for (
int p = 0; p < this->number_of_phases_; ++p) {
1180 ws.surface_rates[p] = -ws.well_potentials[p];
1186 case Well::ProducerCMode::GRUP:
1188 assert(well.isAvailableForGroupControl());
1189 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1190 const double efficiencyFactor = well.getEfficiencyFactor();
1191 double scale = this->getGroupProductionTargetRate(group,
1201 for (
int p = 0; p<np; ++p) {
1202 ws.surface_rates[p] *= scale;
1204 ws.trivial_target =
false;
1206 ws.trivial_target =
true;
1210 case Well::ProducerCMode::CMODE_UNDEFINED:
1211 case Well::ProducerCMode::NONE:
1213 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1220 ws.bhp = controls.bhp_limit;
1225 template<
typename TypeTag>
1227 WellInterface<TypeTag>::
1228 initialWellRateFractions(
const Simulator& ebosSimulator,
const WellState& well_state)
const
1230 const int np = this->number_of_phases_;
1231 std::vector<double> scaling_factor(np);
1232 const auto& ws = well_state.well(this->index_of_well_);
1234 double total_potentials = 0.0;
1235 for (
int p = 0; p<np; ++p) {
1236 total_potentials += ws.well_potentials[p];
1238 if (total_potentials > 0) {
1239 for (
int p = 0; p<np; ++p) {
1240 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1242 return scaling_factor;
1246 double total_tw = 0;
1247 const int nperf = this->number_of_perforations_;
1248 for (
int perf = 0; perf < nperf; ++perf) {
1249 total_tw += this->well_index_[perf];
1251 for (
int perf = 0; perf < nperf; ++perf) {
1252 const int cell_idx = this->well_cells_[perf];
1253 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1254 const auto& fs = intQuants.fluidState();
1255 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1256 double total_mobility = 0.0;
1257 for (
int p = 0; p < np; ++p) {
1258 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1259 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1261 for (
int p = 0; p < np; ++p) {
1262 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1263 scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1266 return scaling_factor;
1271 template <
typename TypeTag>
1280 auto&
ws = well_state.well(this->index_of_well_);
1283 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1300 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1301 ws.surface_rates[
p] =
well_q_s[this->flowPhaseToEbosCompIdx(
p)];
1310 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1312 const int comp_idx = this->flowPhaseToEbosCompIdx(
p);
1313 double&
rate =
ws.surface_rates[
p];
1320 template<
typename TypeTag>
1321 typename WellInterface<TypeTag>::Eval
1324 if constexpr (Indices::oilEnabled) {
1325 return fs.pressure(FluidSystem::oilPhaseIdx);
1326 }
else if constexpr (Indices::gasEnabled) {
1327 return fs.pressure(FluidSystem::gasPhaseIdx);
1329 return fs.pressure(FluidSystem::waterPhaseIdx);
1333 template <
typename TypeTag>
1334 template<
class Value,
class Callback>
1336 WellInterface<TypeTag>::
1337 getMobility(
const Simulator& ebosSimulator,
1339 std::vector<Value>& mob,
1340 Callback& extendEval,
1341 [[maybe_unused]] DeferredLogger& deferred_logger)
const
1343 auto relpermArray = []()
1345 if constexpr (std::is_same_v<Value, Scalar>) {
1346 return std::array<Scalar,3>{};
1348 return std::array<Eval,3>{};
1351 const int cell_idx = this->well_cells_[perf];
1352 assert (
int(mob.size()) == this->num_components_);
1353 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1354 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1358 const int satid = this->saturation_table_number_[perf] - 1;
1359 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1360 if (satid == satid_elem) {
1361 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1362 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1366 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1367 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1369 if constexpr (has_solvent) {
1370 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1373 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1374 auto relativePerms = relpermArray();
1375 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1378 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1381 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1382 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1386 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1387 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1391 if constexpr (has_solvent) {
1392 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1396 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1397 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1398 const auto& connections = this->well_ecl_.getConnections();
1399 const auto& connection = connections[perf_ecl_index];
1400 if (connection.filterCakeActive()) {
1401 for (
auto& val : mob) {
1402 val *= this->inj_fc_multiplier_[perf];
1409 template<
typename TypeTag>
1411 WellInterface<TypeTag>::
1412 updateWellStateWithTHPTargetProd(
const Simulator& ebos_simulator,
1413 WellState& well_state,
1414 DeferredLogger& deferred_logger)
const
1416 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1418 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1419 ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
1420 if (bhp_at_thp_limit) {
1421 std::vector<double> rates(this->number_of_phases_, 0.0);
1422 if (thp_update_iterations) {
1423 computeWellRatesWithBhpIterations(ebos_simulator, *bhp_at_thp_limit,
1424 rates, deferred_logger);
1426 computeWellRatesWithBhp(ebos_simulator, *bhp_at_thp_limit,
1427 rates, deferred_logger);
1429 auto& ws = well_state.well(this->name());
1430 ws.surface_rates = rates;
1431 ws.bhp = *bhp_at_thp_limit;
1432 ws.thp = this->getTHPConstraint(summary_state);
1439 template <
typename TypeTag>
1441 WellInterface<TypeTag>::
1442 computeConnLevelProdInd(
const FluidState& fs,
1443 const std::function<
double(
const double)>& connPICalc,
1444 const std::vector<Scalar>& mobility,
1445 double* connPI)
const
1448 const int np = this->number_of_phases_;
1449 for (
int p = 0; p < np; ++p) {
1452 const auto connMob =
1453 mobility[this->flowPhaseToEbosCompIdx(p)]
1454 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
1456 connPI[p] = connPICalc(connMob);
1459 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1460 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1462 const auto io = pu.phase_pos[Oil];
1463 const auto ig = pu.phase_pos[Gas];
1465 const auto vapoil = connPI[ig] * fs.Rv().value();
1466 const auto disgas = connPI[io] * fs.Rs().value();
1468 connPI[io] += vapoil;
1469 connPI[ig] += disgas;
1474 template <
typename TypeTag>
1476 WellInterface<TypeTag>::
1477 computeConnLevelInjInd(
const FluidState& fs,
1478 const Phase preferred_phase,
1479 const std::function<
double(
const double)>& connIICalc,
1480 const std::vector<Scalar>& mobility,
1482 DeferredLogger& deferred_logger)
const
1488 if (preferred_phase == Phase::GAS) {
1489 phase_pos = pu.phase_pos[Gas];
1491 else if (preferred_phase == Phase::OIL) {
1492 phase_pos = pu.phase_pos[Oil];
1494 else if (preferred_phase == Phase::WATER) {
1495 phase_pos = pu.phase_pos[Water];
1498 OPM_DEFLOG_THROW(NotImplemented,
1499 fmt::format(
"Unsupported Injector Type ({}) "
1500 "for well {} during connection I.I. calculation",
1501 static_cast<int>(preferred_phase), this->name()),
1505 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1506 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
Definition AquiferInterface.hpp:35
Definition DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition ParallelWellInfo.hpp:184
Definition WellInterfaceIndices.hpp:35
Definition WellInterface.hpp:74
void updateWellStateRates(const Simulator &ebosSimulator, WellState &well_state, DeferredLogger &deferred_logger) const
Modify the well_state's rates if there is only one nonzero rate.
Definition WellInterface_impl.hpp:1274
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:43
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
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 BlackoilModelParametersEbos.hpp:451
Definition BlackoilPhases.hpp:46