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 const auto& schedule = simulator.vanguard().schedule();
357 const auto report_step = simulator.episodeIndex();
358 const auto& glo = schedule.glo(report_step);
360 gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
364 WellTestState welltest_state_temp;
366 bool testWell =
true;
371 const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
372 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
374 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
375 deferred_logger.debug(msg);
380 updateWellOperability(simulator, well_state_copy, deferred_logger);
381 if ( !this->isOperableAndSolvable() ) {
382 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
383 deferred_logger.debug(msg);
386 std::vector<double> potentials;
388 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
389 }
catch (
const std::exception& e) {
390 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during testing for re-opening: ") + e.what();
391 deferred_logger.info(msg);
394 const int np = well_state_copy.numPhases();
395 for (
int p = 0; p < np; ++p) {
396 ws.well_potentials[p] = std::max(0.0, potentials[p]);
398 this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time,
false, welltest_state_temp, deferred_logger);
399 this->closeCompletions(welltest_state_temp);
405 if ( welltest_state_temp.num_closed_wells() > 0 ||
406 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
412 if (!welltest_state_temp.well_is_closed(this->name())) {
413 well_test_state.open_well(this->name());
415 std::string msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
416 deferred_logger.info(msg);
419 for (
auto& completion : this->well_ecl_.getCompletions()) {
420 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
421 well_test_state.open_completion(this->name(), completion.first);
425 well_state = well_state_copy;
432 template<
typename TypeTag>
434 WellInterface<TypeTag>::
435 iterateWellEquations(
const Simulator& ebosSimulator,
437 WellState& well_state,
438 const GroupState& group_state,
439 DeferredLogger& deferred_logger)
441 const auto& summary_state = ebosSimulator.vanguard().summaryState();
442 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
443 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
444 bool converged =
false;
447 if (!this->param_.local_well_solver_control_switching_){
448 converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
450 converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
453 }
catch (NumericalProblem& e ) {
454 const std::string msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
455 deferred_logger.warning(
"INNER_ITERATION_FAILED", msg);
462 template<
typename TypeTag>
464 WellInterface<TypeTag>::
465 solveWellForTesting(
const Simulator& ebosSimulator, WellState& well_state,
const GroupState& group_state,
466 DeferredLogger& deferred_logger)
469 const WellState well_state0 = well_state;
470 const double dt = ebosSimulator.timeStepSize();
471 const auto& summary_state = ebosSimulator.vanguard().summaryState();
472 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
475 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
476 converged = gliftBeginTimeStepWellTestIterateWellEquations(
477 ebosSimulator, dt, well_state, group_state, deferred_logger);
480 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
481 converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
484 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
487 const int max_iter = param_.max_welleq_iter_;
488 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
489 + std::to_string(max_iter) +
" iterations");
490 well_state = well_state0;
495 template<
typename TypeTag>
497 WellInterface<TypeTag>::
498 solveWellEquation(
const Simulator& ebosSimulator,
499 WellState& well_state,
500 const GroupState& group_state,
501 DeferredLogger& deferred_logger)
503 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
507 const WellState well_state0 = well_state;
508 const double dt = ebosSimulator.timeStepSize();
509 bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
519 auto& ws = well_state.well(this->indexOfWell());
520 bool thp_control =
false;
521 if (this->well_ecl_.isInjector()) {
522 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
524 ws.injection_cmode = Well::InjectorCMode::BHP;
525 this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
528 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
530 ws.production_cmode = Well::ProducerCMode::BHP;
531 this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
535 const std::string msg = std::string(
"The newly opened well ") + this->name()
536 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
537 deferred_logger.debug(msg);
538 converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
543 const int max_iter = param_.max_welleq_iter_;
544 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
545 + std::to_string(max_iter) +
" iterations");
546 well_state = well_state0;
552 template <
typename TypeTag>
554 WellInterface<TypeTag>::
555 assembleWellEq(
const Simulator& ebosSimulator,
557 WellState& well_state,
558 const GroupState& group_state,
559 DeferredLogger& deferred_logger)
562 prepareWellBeforeAssembling(ebosSimulator, dt, well_state, group_state, deferred_logger);
564 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, group_state, deferred_logger);
569 template <
typename TypeTag>
571 WellInterface<TypeTag>::
572 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
574 WellState& well_state,
575 const GroupState& group_state,
576 DeferredLogger& deferred_logger)
578 const auto& summary_state = ebosSimulator.vanguard().summaryState();
579 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
580 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
583 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
588 template<
typename TypeTag>
590 WellInterface<TypeTag>::
591 prepareWellBeforeAssembling(
const Simulator& ebosSimulator,
593 WellState& well_state,
594 const GroupState& group_state,
595 DeferredLogger& deferred_logger)
597 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
599 if (param_.check_well_operability_iter_)
600 checkWellOperability(ebosSimulator, well_state, deferred_logger);
603 const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
604 if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
605 this->operability_status_.solvable =
true;
606 bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
610 if (param_.shut_unsolvable_wells_)
611 this->operability_status_.solvable =
false;
614 if (this->operability_status_.has_negative_potentials) {
615 auto well_state_copy = well_state;
616 std::vector<double> potentials;
618 computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
619 }
catch (
const std::exception& e) {
620 const std::string msg = std::string(
"well ") + this->name() + std::string(
": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
621 deferred_logger.info(msg);
622 this->operability_status_.has_negative_potentials =
true;
624 auto& ws = well_state.well(this->indexOfWell());
625 const int np = well_state.numPhases();
626 for (
int p = 0; p < np; ++p) {
627 ws.well_potentials[p] = std::max(0.0, potentials[p]);
630 this->changed_to_open_this_step_ =
false;
631 const bool well_operable = this->operability_status_.isOperableAndSolvable();
633 if (!well_operable && old_well_operable) {
634 if (this->param_.local_well_solver_control_switching_) {
635 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
637 changed_to_stopped_this_step_ =
true;
641 if (this->well_ecl_.getAutomaticShutIn()) {
642 deferred_logger.info(
" well " + this->name() +
" gets SHUT during iteration ");
644 if (!this->wellIsStopped()) {
645 deferred_logger.info(
" well " + this->name() +
" gets STOPPED during iteration ");
647 changed_to_stopped_this_step_ =
true;
651 }
else if (well_operable && !old_well_operable) {
652 deferred_logger.info(
" well " + this->name() +
" gets REVIVED during iteration ");
654 changed_to_stopped_this_step_ =
false;
655 this->changed_to_open_this_step_ =
true;
659 template<
typename TypeTag>
661 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
663 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
666 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
667 if (this->cells()[perfIdx] == cellIdx) {
668 for (
int i = 0; i < RateVector::dimension; ++i) {
669 rates[i] += connectionRates_[perfIdx][i];
675 template<
typename TypeTag>
676 typename WellInterface<TypeTag>::Scalar
677 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const {
678 for (
int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
679 if (this->cells()[perfIdx] == cellIdx) {
680 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
681 return connectionRates_[perfIdx][activeCompIdx].value();
685 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
686 +
" does not perforate cell " + std::to_string(cellIdx));
693 template<
typename TypeTag>
695 WellInterface<TypeTag>::
696 checkWellOperability(
const Simulator& ebos_simulator,
697 const WellState& well_state,
698 DeferredLogger& deferred_logger)
701 if (!param_.check_well_operability_) {
705 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
709 updateWellOperability(ebos_simulator, well_state, deferred_logger);
710 if (!this->operability_status_.isOperableAndSolvable()) {
711 this->operability_status_.use_vfpexplicit =
true;
712 deferred_logger.debug(
"EXPLICIT_LOOKUP_VFP",
713 "well not operable, trying with explicit vfp lookup: " + this->name());
714 updateWellOperability(ebos_simulator, well_state, deferred_logger);
718 template<
typename TypeTag>
720 WellInterface<TypeTag>::
721 gliftBeginTimeStepWellTestIterateWellEquations(
722 const Simulator& ebos_simulator,
724 WellState& well_state,
725 const GroupState &group_state,
726 DeferredLogger& deferred_logger)
728 const auto& well_name = this->name();
729 assert(this->wellHasTHPConstraints(ebos_simulator.vanguard().summaryState()));
730 const auto& schedule = ebos_simulator.vanguard().schedule();
731 auto report_step_idx = ebos_simulator.episodeIndex();
732 const auto& glo = schedule.glo(report_step_idx);
733 if(glo.active() && glo.has_well(well_name)) {
734 const auto increment = glo.gaslift_increment();
735 auto alq = well_state.getALQ(well_name);
738 well_state.setALQ(well_name, alq);
740 iterateWellEquations(ebos_simulator, dt, well_state, group_state, deferred_logger)))
749 return iterateWellEquations(ebos_simulator, dt, well_state, group_state, deferred_logger);
753 template<
typename TypeTag>
755 WellInterface<TypeTag>::
756 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& ebos_simulator,
757 WellState& well_state,
758 DeferredLogger& deferred_logger)
760 const auto& summary_state = ebos_simulator.vanguard().summaryState();
761 const auto& well_name = this->name();
762 if (!this->wellHasTHPConstraints(summary_state)) {
763 const std::string msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
764 deferred_logger.info(msg);
767 const auto& schedule = ebos_simulator.vanguard().schedule();
768 const auto report_step_idx = ebos_simulator.episodeIndex();
769 const auto& glo = schedule.glo(report_step_idx);
770 if (!glo.has_well(well_name)) {
771 const std::string msg = fmt::format(
772 "GLIFT WTEST: Well {} : Gas Lift not activated: "
773 "WLIFTOPT is probably missing. Skipping.", well_name);
774 deferred_logger.info(msg);
777 const auto& gl_well = glo.well(well_name);
778 auto& max_alq_optional = gl_well.max_rate();
780 if (max_alq_optional) {
781 max_alq = *max_alq_optional;
784 const auto& well_ecl = this->wellEcl();
785 const auto& controls = well_ecl.productionControls(summary_state);
786 const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
787 const auto& alq_values = table.getALQAxis();
788 max_alq = alq_values.back();
790 well_state.setALQ(well_name, max_alq);
791 const std::string msg = fmt::format(
792 "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
794 deferred_logger.info(msg);
797 template<
typename TypeTag>
799 WellInterface<TypeTag>::
800 updateWellOperability(
const Simulator& ebos_simulator,
801 const WellState& well_state,
802 DeferredLogger& deferred_logger)
804 if (this->param_.local_well_solver_control_switching_) {
805 const bool success = updateWellOperabilityFromWellEq(ebos_simulator, well_state, deferred_logger);
809 deferred_logger.debug(
"Operability check using well equations did not converge for well "
810 + this->name() +
", reverting to classical approach." );
813 this->operability_status_.resetOperability();
815 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
816 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
817 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
818 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
822 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
823 if (check_thp || bhp_controlled) {
824 updateIPR(ebos_simulator, deferred_logger);
825 checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
829 checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
833 template<
typename TypeTag>
835 WellInterface<TypeTag>::
836 updateWellOperabilityFromWellEq(
const Simulator& ebos_simulator,
837 const WellState& well_state,
838 DeferredLogger& deferred_logger)
841 assert(this->param_.local_well_solver_control_switching_);
842 this->operability_status_.resetOperability();
843 WellState well_state_copy = well_state;
844 const auto& group_state = ebos_simulator.problem().wellModel().groupState();
845 const double dt = ebos_simulator.timeStepSize();
847 bool converged = iterateWellEquations(ebos_simulator, dt, well_state_copy, group_state, deferred_logger);
851 template<
typename TypeTag>
853 WellInterface<TypeTag>::
854 updateWellStateWithTarget(
const Simulator& ebos_simulator,
855 const GroupState& group_state,
856 WellState& well_state,
857 DeferredLogger& deferred_logger)
const
861 const auto& well = this->well_ecl_;
862 const int well_index = this->index_of_well_;
863 auto& ws = well_state.well(well_index);
865 const int np = well_state.numPhases();
866 const auto& summaryState = ebos_simulator.vanguard().summaryState();
867 const auto& schedule = ebos_simulator.vanguard().schedule();
869 if (this->wellIsStopped()) {
870 for (
int p = 0; p<np; ++p) {
871 ws.surface_rates[p] = 0;
877 if (this->isInjector() )
879 const auto& controls = well.injectionControls(summaryState);
881 InjectorType injectorType = controls.injector_type;
883 switch (injectorType) {
884 case InjectorType::WATER:
886 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
889 case InjectorType::OIL:
891 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
894 case InjectorType::GAS:
896 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
900 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
903 const auto current = ws.injection_cmode;
906 case Well::InjectorCMode::RATE:
908 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
909 if(this->rsRvInj() > 0) {
910 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
911 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
912 }
else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
913 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
915 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 );
921 case Well::InjectorCMode::RESV:
923 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
924 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, convert_coeff);
925 const double coeff = convert_coeff[phasePos];
926 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
930 case Well::InjectorCMode::THP:
932 auto rates = ws.surface_rates;
933 double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
937 this->getRefDensity(),
940 ws.thp = this->getTHPConstraint(summaryState);
945 double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
946 if (total_rate <= 0.0)
947 ws.surface_rates = ws.well_potentials;
951 case Well::InjectorCMode::BHP:
953 ws.bhp = controls.bhp_limit;
954 double total_rate = 0.0;
955 for (
int p = 0; p<np; ++p) {
956 total_rate += ws.surface_rates[p];
961 if (total_rate <= 0.0)
962 ws.surface_rates = ws.well_potentials;
966 case Well::InjectorCMode::GRUP:
968 assert(well.isAvailableForGroupControl());
969 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
970 const double efficiencyFactor = well.getEfficiencyFactor();
971 std::optional<double> target =
972 this->getGroupInjectionTargetRate(group,
981 ws.surface_rates[phasePos] = *target;
984 case Well::InjectorCMode::CMODE_UNDEFINED:
986 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(), deferred_logger );
996 ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
999 ws.bhp = controls.bhp_limit;
1005 const auto current = ws.production_cmode;
1006 const auto& controls = well.productionControls(summaryState);
1008 case Well::ProducerCMode::ORAT:
1010 double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
1013 if (current_rate > 0.0) {
1014 for (
int p = 0; p<np; ++p) {
1015 ws.surface_rates[p] *= controls.oil_rate/current_rate;
1018 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1019 double control_fraction = fractions[pu.phase_pos[Oil]];
1020 if (control_fraction != 0.0) {
1021 for (
int p = 0; p<np; ++p) {
1022 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1028 case Well::ProducerCMode::WRAT:
1030 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
1033 if (current_rate > 0.0) {
1034 for (
int p = 0; p<np; ++p) {
1035 ws.surface_rates[p] *= controls.water_rate/current_rate;
1038 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1039 double control_fraction = fractions[pu.phase_pos[Water]];
1040 if (control_fraction != 0.0) {
1041 for (
int p = 0; p<np; ++p) {
1042 ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
1048 case Well::ProducerCMode::GRAT:
1050 double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
1053 if (current_rate > 0.0) {
1054 for (
int p = 0; p<np; ++p) {
1055 ws.surface_rates[p] *= controls.gas_rate/current_rate;
1058 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1059 double control_fraction = fractions[pu.phase_pos[Gas]];
1060 if (control_fraction != 0.0) {
1061 for (
int p = 0; p<np; ++p) {
1062 ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
1070 case Well::ProducerCMode::LRAT:
1072 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
1073 - ws.surface_rates[ pu.phase_pos[Oil] ];
1076 if (current_rate > 0.0) {
1077 for (
int p = 0; p<np; ++p) {
1078 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1081 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1082 double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
1083 if (control_fraction != 0.0) {
1084 for (
int p = 0; p<np; ++p) {
1085 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1091 case Well::ProducerCMode::CRAT:
1093 OPM_DEFLOG_THROW(std::runtime_error,
1094 fmt::format(
"CRAT control not supported, well {}", this->name()),
1097 case Well::ProducerCMode::RESV:
1099 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
1100 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1101 double total_res_rate = 0.0;
1102 for (
int p = 0; p<np; ++p) {
1103 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1105 if (controls.prediction_mode) {
1108 if (total_res_rate > 0.0) {
1109 for (
int p = 0; p<np; ++p) {
1110 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1113 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1114 for (
int p = 0; p<np; ++p) {
1115 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1119 std::vector<double> hrates(this->number_of_phases_,0.);
1120 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1121 hrates[pu.phase_pos[Water]] = controls.water_rate;
1123 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1124 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
1126 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1127 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
1129 std::vector<double> hrates_resv(this->number_of_phases_,0.);
1130 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_, hrates, hrates_resv);
1131 double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1134 if (total_res_rate > 0.0) {
1135 for (
int p = 0; p<np; ++p) {
1136 ws.surface_rates[p] *= target/total_res_rate;
1139 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1140 for (
int p = 0; p<np; ++p) {
1141 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1148 case Well::ProducerCMode::BHP:
1150 ws.bhp = controls.bhp_limit;
1151 double total_rate = 0.0;
1152 for (
int p = 0; p<np; ++p) {
1153 total_rate -= ws.surface_rates[p];
1158 if (total_rate <= 0.0){
1159 for (
int p = 0; p<np; ++p) {
1160 ws.surface_rates[p] = -ws.well_potentials[p];
1165 case Well::ProducerCMode::THP:
1167 const bool update_success = updateWellStateWithTHPTargetProd(ebos_simulator, well_state, deferred_logger);
1169 if (!update_success) {
1173 auto rates = ws.surface_rates;
1174 this->adaptRatesForVFP(rates);
1175 const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1176 well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1178 ws.thp = this->getTHPConstraint(summaryState);
1182 const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1183 if (total_rate <= 0.0) {
1184 for (
int p = 0; p < this->number_of_phases_; ++p) {
1185 ws.surface_rates[p] = -ws.well_potentials[p];
1191 case Well::ProducerCMode::GRUP:
1193 assert(well.isAvailableForGroupControl());
1194 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1195 const double efficiencyFactor = well.getEfficiencyFactor();
1196 double scale = this->getGroupProductionTargetRate(group,
1206 for (
int p = 0; p<np; ++p) {
1207 ws.surface_rates[p] *= scale;
1209 ws.trivial_target =
false;
1211 ws.trivial_target =
true;
1215 case Well::ProducerCMode::CMODE_UNDEFINED:
1216 case Well::ProducerCMode::NONE:
1218 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() , deferred_logger);
1225 ws.bhp = controls.bhp_limit;
1230 template<
typename TypeTag>
1232 WellInterface<TypeTag>::
1233 initialWellRateFractions(
const Simulator& ebosSimulator,
const WellState& well_state)
const
1235 const int np = this->number_of_phases_;
1236 std::vector<double> scaling_factor(np);
1237 const auto& ws = well_state.well(this->index_of_well_);
1239 double total_potentials = 0.0;
1240 for (
int p = 0; p<np; ++p) {
1241 total_potentials += ws.well_potentials[p];
1243 if (total_potentials > 0) {
1244 for (
int p = 0; p<np; ++p) {
1245 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1247 return scaling_factor;
1251 double total_tw = 0;
1252 const int nperf = this->number_of_perforations_;
1253 for (
int perf = 0; perf < nperf; ++perf) {
1254 total_tw += this->well_index_[perf];
1256 for (
int perf = 0; perf < nperf; ++perf) {
1257 const int cell_idx = this->well_cells_[perf];
1258 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1259 const auto& fs = intQuants.fluidState();
1260 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1261 double total_mobility = 0.0;
1262 for (
int p = 0; p < np; ++p) {
1263 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1264 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1266 for (
int p = 0; p < np; ++p) {
1267 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1268 scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1271 return scaling_factor;
1276 template <
typename TypeTag>
1285 auto&
ws = well_state.well(this->index_of_well_);
1288 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1305 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1306 ws.surface_rates[
p] =
well_q_s[this->flowPhaseToEbosCompIdx(
p)];
1315 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1317 const int comp_idx = this->flowPhaseToEbosCompIdx(
p);
1318 double&
rate =
ws.surface_rates[
p];
1325 template<
typename TypeTag>
1326 typename WellInterface<TypeTag>::Eval
1329 if constexpr (Indices::oilEnabled) {
1330 return fs.pressure(FluidSystem::oilPhaseIdx);
1331 }
else if constexpr (Indices::gasEnabled) {
1332 return fs.pressure(FluidSystem::gasPhaseIdx);
1334 return fs.pressure(FluidSystem::waterPhaseIdx);
1338 template <
typename TypeTag>
1339 template<
class Value,
class Callback>
1341 WellInterface<TypeTag>::
1342 getMobility(
const Simulator& ebosSimulator,
1344 std::vector<Value>& mob,
1345 Callback& extendEval,
1346 [[maybe_unused]] DeferredLogger& deferred_logger)
const
1348 auto relpermArray = []()
1350 if constexpr (std::is_same_v<Value, Scalar>) {
1351 return std::array<Scalar,3>{};
1353 return std::array<Eval,3>{};
1356 const int cell_idx = this->well_cells_[perf];
1357 assert (
int(mob.size()) == this->num_components_);
1358 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1359 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1363 const int satid = this->saturation_table_number_[perf] - 1;
1364 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1365 if (satid == satid_elem) {
1366 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1367 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1371 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1372 mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
1374 if constexpr (has_solvent) {
1375 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1378 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1379 auto relativePerms = relpermArray();
1380 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1383 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1386 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1387 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1391 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1392 mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1396 if constexpr (has_solvent) {
1397 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent", deferred_logger);
1401 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1402 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1403 const auto& connections = this->well_ecl_.getConnections();
1404 const auto& connection = connections[perf_ecl_index];
1405 if (connection.filterCakeActive()) {
1406 for (
auto& val : mob) {
1407 val *= this->inj_fc_multiplier_[perf];
1414 template<
typename TypeTag>
1416 WellInterface<TypeTag>::
1417 updateWellStateWithTHPTargetProd(
const Simulator& ebos_simulator,
1418 WellState& well_state,
1419 DeferredLogger& deferred_logger)
const
1421 const auto& summary_state = ebos_simulator.vanguard().summaryState();
1423 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1424 ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
1425 if (bhp_at_thp_limit) {
1426 std::vector<double> rates(this->number_of_phases_, 0.0);
1427 if (thp_update_iterations) {
1428 computeWellRatesWithBhpIterations(ebos_simulator, *bhp_at_thp_limit,
1429 rates, deferred_logger);
1431 computeWellRatesWithBhp(ebos_simulator, *bhp_at_thp_limit,
1432 rates, deferred_logger);
1434 auto& ws = well_state.well(this->name());
1435 ws.surface_rates = rates;
1436 ws.bhp = *bhp_at_thp_limit;
1437 ws.thp = this->getTHPConstraint(summary_state);
1444 template <
typename TypeTag>
1446 WellInterface<TypeTag>::
1447 computeConnLevelProdInd(
const FluidState& fs,
1448 const std::function<
double(
const double)>& connPICalc,
1449 const std::vector<Scalar>& mobility,
1450 double* connPI)
const
1453 const int np = this->number_of_phases_;
1454 for (
int p = 0; p < np; ++p) {
1457 const auto connMob =
1458 mobility[this->flowPhaseToEbosCompIdx(p)]
1459 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
1461 connPI[p] = connPICalc(connMob);
1464 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1465 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1467 const auto io = pu.phase_pos[Oil];
1468 const auto ig = pu.phase_pos[Gas];
1470 const auto vapoil = connPI[ig] * fs.Rv().value();
1471 const auto disgas = connPI[io] * fs.Rs().value();
1473 connPI[io] += vapoil;
1474 connPI[ig] += disgas;
1479 template <
typename TypeTag>
1481 WellInterface<TypeTag>::
1482 computeConnLevelInjInd(
const FluidState& fs,
1483 const Phase preferred_phase,
1484 const std::function<
double(
const double)>& connIICalc,
1485 const std::vector<Scalar>& mobility,
1487 DeferredLogger& deferred_logger)
const
1493 if (preferred_phase == Phase::GAS) {
1494 phase_pos = pu.phase_pos[Gas];
1496 else if (preferred_phase == Phase::OIL) {
1497 phase_pos = pu.phase_pos[Oil];
1499 else if (preferred_phase == Phase::WATER) {
1500 phase_pos = pu.phase_pos[Water];
1503 OPM_DEFLOG_THROW(NotImplemented,
1504 fmt::format(
"Unsupported Injector Type ({}) "
1505 "for well {} during connection I.I. calculation",
1506 static_cast<int>(preferred_phase), this->name()),
1510 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1511 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:1279
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