22#ifndef OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
23#define OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
26#ifndef OPM_WELLINTERFACE_HEADER_INCLUDED
28#include <opm/simulators/wells/WellInterface.hpp>
31#include <opm/common/Exceptions.hpp>
33#include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
34#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
36#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
38#include <opm/simulators/wells/GroupState.hpp>
39#include <opm/simulators/wells/TargetCalculator.hpp>
40#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
41#include <opm/simulators/wells/WellHelpers.hpp>
43#include <dune/common/version.hh>
50#include <fmt/format.h>
56 template<
typename TypeTag>
61 const ModelParameters& param,
63 const int pvtRegionIdx,
79 connectionRates_.resize(this->number_of_local_perforations_);
81 if constexpr (has_solvent || has_zFraction) {
82 if (well.isInjector()) {
85 this->wsolvent_ = this->well_ecl_.getSolventFraction();
92 template<
typename TypeTag>
96 const std::vector<Scalar>& ,
98 const std::vector<Scalar>&
B_avg,
110 template<
typename TypeTag>
111 typename WellInterface<TypeTag>::Scalar
112 WellInterface<TypeTag>::
115 if constexpr (has_polymer) {
116 return this->wpolymer_();
126 template<
typename TypeTag>
127 typename WellInterface<TypeTag>::Scalar
128 WellInterface<TypeTag>::
131 if constexpr (has_foam) {
132 return this->wfoam_();
140 template<
typename TypeTag>
141 typename WellInterface<TypeTag>::Scalar
142 WellInterface<TypeTag>::
145 if constexpr (has_brine) {
146 return this->wsalt_();
152 template<
typename TypeTag>
153 typename WellInterface<TypeTag>::Scalar
154 WellInterface<TypeTag>::
157 if constexpr (has_micp) {
158 return this->wmicrobes_();
164 template<
typename TypeTag>
165 typename WellInterface<TypeTag>::Scalar
166 WellInterface<TypeTag>::
169 if constexpr (has_micp) {
170 return this->woxygen_();
182 template<
typename TypeTag>
183 typename WellInterface<TypeTag>::Scalar
184 WellInterface<TypeTag>::
187 if constexpr (has_micp) {
188 return this->wurea_();
194 template<
typename TypeTag>
196 WellInterface<TypeTag>::
197 updateWellControl(
const Simulator& simulator,
198 const IndividualOrGroup
iog,
208 const auto& summaryState = simulator.vanguard().summaryState();
209 const auto& schedule = simulator.vanguard().schedule();
210 const auto& well = this->well_ecl_;
211 auto&
ws = well_state.well(this->index_of_well_);
213 if (well.isInjector()) {
218 bool oscillating = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) >= this->param_.max_number_of_well_switches_;
219 const int episodeIdx = simulator.episodeIndex();
220 const int iterationIdx = simulator.model().newtonMethod().numIterations();
224 bool output = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) == this->param_.max_number_of_well_switches_;
226 std::ostringstream
ss;
227 ss <<
" The control mode for well " << this->name()
228 <<
" is oscillating\n"
229 <<
" We don't allow for more than "
230 << this->param_.max_number_of_well_switches_
231 <<
" switches. The control is kept at " <<
from;
234 this->well_control_log_.push_back(
from);
239 if (
iog == IndividualOrGroup::Individual) {
241 }
else if (
iog == IndividualOrGroup::Group) {
247 Parallel::Communication
cc = simulator.vanguard().grid().comm();
251 if (well.isInjector()) {
256 std::ostringstream
ss;
257 ss <<
" Switching control mode for well " << this->name()
261 ss <<
" on rank " <<
cc.rank();
265 this->well_control_log_.push_back(
from);
266 updateWellStateWithTarget(simulator, group_state, well_state,
deferred_logger);
273 template<
typename TypeTag>
275 WellInterface<TypeTag>::
276 updateWellControlAndStatusLocalIteration(
const Simulator& simulator,
287 const auto&
summary_state = simulator.vanguard().summaryState();
288 const auto& schedule = simulator.vanguard().schedule();
289 auto&
ws = well_state.well(this->index_of_well_);
291 if (this->isInjector()) {
296 const bool oscillating = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) >= this->param_.max_number_of_well_switches_;
298 if (
oscillating || this->wellUnderZeroRateTarget(simulator, well_state,
deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
302 const Scalar
sgn = this->isInjector() ? 1.0 : -1.0;
303 if (!this->wellIsStopped()){
316 bool isGroupControl =
ws.production_cmode == Well::ProducerCMode::GRUP ||
ws.injection_cmode == Well::InjectorCMode::GRUP;
317 if (! (
isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
320 if (
hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
325 const bool thp_controlled = this->isInjector() ?
ws.injection_cmode == Well::InjectorCMode::THP :
326 ws.production_cmode == Well::ProducerCMode::THP;
329 updateWellStateWithTarget(simulator, group_state, well_state,
deferred_logger);
340 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
345 std::vector<Scalar> rates(this->num_components_);
346 if (this->isInjector()){
347 const Scalar
bhp_thp = WellBhpThpCalculator(*this).
348 calculateBhpFromThp(well_state, rates,
351 this->getRefDensity(),
357 const Scalar
bhp_min = WellBhpThpCalculator(*this).
358 calculateMinimumBhpFromThp(well_state,
361 this->getRefDensity());
370 well_state.well(this->index_of_well_).thp = this->getTHPConstraint(
summary_state);
381 template<
typename TypeTag>
383 WellInterface<TypeTag>::
384 wellTesting(
const Simulator& simulator,
400 const auto&
summary_state = simulator.vanguard().summaryState();
403 ws.production_cmode = Well::ProducerCMode::THP;
406 ws.production_cmode = Well::ProducerCMode::BHP;
413 if (this->isProducer()) {
414 const auto& schedule = simulator.vanguard().schedule();
415 const auto report_step = simulator.episodeIndex();
416 const auto&
glo = schedule.glo(report_step);
418 gliftBeginTimeStepWellTestUpdateALQ(simulator,
437 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
444 if ( !this->isOperableAndSolvable() ) {
445 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
452 }
catch (
const std::exception&
e) {
453 const std::string
msg = fmt::format(
"well {}: computeWellPotentials() "
454 "failed during testing for re-opening: ",
455 this->name(),
e.what());
460 for (
int p = 0;
p <
np; ++
p) {
484 well_test_state.open_well(this->name());
486 std::string
msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
492 well_test_state.open_completion(this->name(),
completion.first);
497 open_times.try_emplace(this->name(), well_test_state.lastTestTime(
this->name()));
504 template<
typename TypeTag>
506 WellInterface<TypeTag>::
507 iterateWellEquations(
const Simulator& simulator,
514 const auto&
summary_state = simulator.vanguard().summaryState();
515 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(
summary_state) : Well::InjectionControls(0);
516 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(
summary_state) : Well::ProductionControls(0);
517 bool converged =
false;
520 if (!this->param_.local_well_solver_control_switching_){
523 if (this->param_.use_implicit_ipr_ &&
this->well_ecl_.isProducer() &&
this->wellHasTHPConstraints(
summary_state) && (
this->well_ecl_.getStatus() == WellStatus::OPEN)) {
531 const std::string
msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
538 template<
typename TypeTag>
540 WellInterface<TypeTag>::
541 solveWellWithTHPConstraint(
const Simulator& simulator,
550 const auto&
summary_state = simulator.vanguard().summaryState();
552 bool converged =
true;
553 auto&
ws = well_state.well(this->index_of_well_);
555 if (this->wellIsStopped()) {
560 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
569 const Scalar bhp = std::max(
bhp_target.value(),
579 const bool isThp =
ws.production_cmode == Well::ProducerCMode::THP;
582 auto rates = well_state.well(this->index_of_well_).surface_rates;
583 this->adaptRatesForVFP(rates);
588 this->operability_status_.use_vfpexplicit =
true;
591 const Scalar reltol = 1
e-3;
594 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
606 this->operability_status_.use_vfpexplicit =
true;
617 const Scalar bhp = std::max(
bhp_target.value(),
621 converged = this->iterateWellEqWithSwitching(simulator,
dt,
631 this->operability_status_.can_obtain_bhp_with_thp_limit =
is_operable;
632 this->operability_status_.obey_thp_limit_under_bhp_limit =
is_operable;
636 template<
typename TypeTag>
637 std::optional<typename WellInterface<TypeTag>::Scalar>
638 WellInterface<TypeTag>::
639 estimateOperableBhp(
const Simulator& simulator,
648 Scalar
bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state,
this->well_ecl_,
summary_state,
this->getRefDensity());
651 if (!converged || this->wellIsStopped()) {
655 auto rates = well_state.well(this->index_of_well_).surface_rates;
656 this->adaptRatesForVFP(rates);
657 return WellBhpThpCalculator(*this).estimateStableBhp(well_state,
this->well_ecl_, rates,
this->getRefDensity(),
summary_state);
660 template<
typename TypeTag>
662 WellInterface<TypeTag>::
663 solveWellWithBhp(
const Simulator& simulator,
674 auto&
ws = well_state.well(this->index_of_well_);
677 if (this->isInjector()) {
681 ws.injection_cmode = Well::InjectorCMode::BHP;
686 ws.production_cmode = Well::ProducerCMode::BHP;
697 template<
typename TypeTag>
699 WellInterface<TypeTag>::
700 solveWellWithZeroRate(
const Simulator& simulator,
718 template<
typename TypeTag>
720 WellInterface<TypeTag>::
721 solveWellForTesting(
const Simulator& simulator,
727 const double dt = simulator.timeStepSize();
728 const bool converged = iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
730 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
733 const int max_iter = this->param_.max_welleq_iter_;
734 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
735 + std::to_string(
max_iter) +
" iterations");
740 template<
typename TypeTag>
742 WellInterface<TypeTag>::
743 solveWellEquation(
const Simulator& simulator,
749 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
754 const double dt = simulator.timeStepSize();
755 bool converged = iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
765 auto&
ws = well_state.well(this->indexOfWell());
767 if (this->well_ecl_.isInjector()) {
768 thp_control =
ws.injection_cmode == Well::InjectorCMode::THP;
770 ws.injection_cmode = Well::InjectorCMode::BHP;
774 thp_control =
ws.production_cmode == Well::ProducerCMode::THP;
776 ws.production_cmode = Well::ProducerCMode::BHP;
781 const std::string
msg = std::string(
"The newly opened well ") + this->name()
782 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
784 converged = this->iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
789 const int max_iter = this->param_.max_welleq_iter_;
790 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
791 + std::to_string(
max_iter) +
" iterations");
798 template <
typename TypeTag>
800 WellInterface<TypeTag>::
801 assembleWellEq(
const Simulator& simulator,
808 prepareWellBeforeAssembling(simulator,
dt, well_state, group_state,
deferred_logger);
809 assembleWellEqWithoutIteration(simulator,
dt, well_state, group_state,
deferred_logger);
814 template <
typename TypeTag>
816 WellInterface<TypeTag>::
817 assembleWellEqWithoutIteration(
const Simulator& simulator,
824 const auto&
summary_state = simulator.vanguard().summaryState();
825 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(
summary_state) : Well::InjectionControls(0);
826 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(
summary_state) : Well::ProductionControls(0);
834 template<
typename TypeTag>
836 WellInterface<TypeTag>::
837 prepareWellBeforeAssembling(
const Simulator& simulator,
844 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
846 if (this->param_.check_well_operability_iter_)
850 const int iteration_idx = simulator.model().newtonMethod().numIterations();
852 const auto&
ws = well_state.well(this->indexOfWell());
857 std::any_of(
ws.surface_rates.begin(),
858 ws.surface_rates.begin() + well_state.numPhases(),
859 [](Scalar rate) { return rate != Scalar(0.0); });
861 this->operability_status_.solvable =
true;
862 bool converged = this->iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
870 this->operability_status_.resetOperability();
872 deferred_logger.debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
877 if (this->isInjector()) {
884 deferred_logger.debug(
" " + this->name() +
" switched from " +
from +
" to " +
to +
" during local solve");
889 if (this->param_.shut_unsolvable_wells_)
890 this->operability_status_.solvable =
false;
893 if (this->operability_status_.has_negative_potentials) {
898 }
catch (
const std::exception&
e) {
899 const std::string
msg = fmt::format(
"well {}: computeWellPotentials() failed "
900 "during attempt to recompute potentials for well: ",
901 this->name(),
e.what());
903 this->operability_status_.has_negative_potentials =
true;
905 auto&
ws = well_state.well(this->indexOfWell());
906 const int np = well_state.numPhases();
907 for (
int p = 0;
p <
np; ++
p) {
911 this->changed_to_open_this_step_ =
false;
912 const bool well_operable = this->operability_status_.isOperableAndSolvable();
915 deferred_logger.debug(
" well " + this->name() +
" gets STOPPED during iteration ");
917 changed_to_stopped_this_step_ =
true;
919 deferred_logger.debug(
" well " + this->name() +
" gets REVIVED during iteration ");
921 changed_to_stopped_this_step_ =
false;
922 this->changed_to_open_this_step_ =
true;
926 template<
typename TypeTag>
928 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
930 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
935 for (
int i = 0; i < RateVector::dimension; ++i) {
936 rates[i] += connectionRates_[
perfIdx][i];
942 template<
typename TypeTag>
943 typename WellInterface<TypeTag>::Scalar
944 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const
948 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
953 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
954 +
" does not perforate cell " + std::to_string(
cellIdx));
961 template<
typename TypeTag>
963 WellInterface<TypeTag>::
964 checkWellOperability(
const Simulator& simulator,
969 if (!this->param_.check_well_operability_) {
973 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
978 if (!this->operability_status_.isOperableAndSolvable()) {
979 this->operability_status_.use_vfpexplicit =
true;
981 "well not operable, trying with explicit vfp lookup: " + this->name());
986 template<
typename TypeTag>
988 WellInterface<TypeTag>::
989 gliftBeginTimeStepWellTestIterateWellEquations(
const Simulator& simulator,
996 const auto& well_name = this->name();
997 assert(this->wellHasTHPConstraints(simulator.vanguard().summaryState()));
998 const auto& schedule = simulator.vanguard().schedule();
1001 if(
glo.active() &&
glo.has_well(well_name)) {
1002 const auto increment =
glo.gaslift_increment();
1003 auto alq = well_state.getALQ(well_name);
1006 well_state.setALQ(well_name, alq);
1017 return iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
1021 template<
typename TypeTag>
1023 WellInterface<TypeTag>::
1024 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& simulator,
1032 const auto&
summary_state = simulator.vanguard().summaryState();
1033 const auto& well_name = this->name();
1035 const std::string
msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1039 const auto& schedule = simulator.vanguard().schedule();
1042 if (!
glo.has_well(well_name)) {
1043 const std::string
msg = fmt::format(
1044 "GLIFT WTEST: Well {} : Gas lift not activated: "
1045 "WLIFTOPT is probably missing. Skipping.", well_name);
1052 std::unique_ptr<GasLiftSingleWell>
glift =
1063 well_state.setALQ(well_name,
wtest_alq);
1065 "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1071 "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1073 unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.getALQ(well_name)));
1078 "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1085 template<
typename TypeTag>
1087 WellInterface<TypeTag>::
1088 updateWellOperability(
const Simulator& simulator,
1093 if (this->param_.local_well_solver_control_switching_) {
1094 const bool success = updateWellOperabilityFromWellEq(simulator, well_state,
deferred_logger);
1098 deferred_logger.debug(
"Operability check using well equations did not converge for well "
1099 + this->name() +
", reverting to classical approach." );
1102 this->operability_status_.resetOperability();
1104 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1105 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1106 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1107 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1114 checkOperabilityUnderBHPLimit(well_state, simulator,
deferred_logger);
1118 checkOperabilityUnderTHPLimit(simulator, well_state,
deferred_logger);
1122 template<
typename TypeTag>
1124 WellInterface<TypeTag>::
1125 updateWellOperabilityFromWellEq(
const Simulator& simulator,
1131 assert(this->param_.local_well_solver_control_switching_);
1132 this->operability_status_.resetOperability();
1134 const auto& group_state = simulator.problem().wellModel().groupState();
1135 const double dt = simulator.timeStepSize();
1141 template<
typename TypeTag>
1143 WellInterface<TypeTag>::
1144 updateWellStateWithTarget(
const Simulator& simulator,
1151 const auto& well = this->well_ecl_;
1152 const int well_index = this->index_of_well_;
1153 auto&
ws = well_state.well(well_index);
1155 const int np = well_state.numPhases();
1156 const auto& summaryState = simulator.vanguard().summaryState();
1157 const auto& schedule = simulator.vanguard().schedule();
1159 if (this->wellIsStopped()) {
1160 for (
int p = 0;
p<
np; ++
p) {
1161 ws.surface_rates[
p] = 0;
1167 if (this->isInjector() )
1169 const auto&
controls = well.injectionControls(summaryState);
1174 case InjectorType::WATER:
1176 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
1179 case InjectorType::OIL:
1181 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
1184 case InjectorType::GAS:
1186 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
1190 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(),
deferred_logger );
1193 const auto current =
ws.injection_cmode;
1196 case Well::InjectorCMode::RATE:
1199 if(this->rsRvInj() > 0) {
1200 if (
injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1201 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] =
controls.surface_rate * this->rsRvInj();
1202 }
else if (
injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1203 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] =
controls.surface_rate * this->rsRvInj();
1205 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 );
1211 case Well::InjectorCMode::RESV:
1213 std::vector<Scalar>
convert_coeff(this->number_of_phases_, 1.0);
1214 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_,
convert_coeff);
1220 case Well::InjectorCMode::THP:
1222 auto rates =
ws.surface_rates;
1223 Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
1227 this->getRefDensity(),
1230 ws.thp = this->getTHPConstraint(summaryState);
1235 Scalar
total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1237 ws.surface_rates =
ws.well_potentials;
1241 case Well::InjectorCMode::BHP:
1245 for (
int p = 0;
p<
np; ++
p) {
1252 ws.surface_rates =
ws.well_potentials;
1256 case Well::InjectorCMode::GRUP:
1258 assert(well.isAvailableForGroupControl());
1259 const auto& group = schedule.getGroup(well.groupName(),
this->currentStep());
1261 well_state[well.name()].efficiency_scaling_factor;
1262 std::optional<Scalar>
target =
1263 this->getGroupInjectionTargetRate(group,
1275 case Well::InjectorCMode::CMODE_UNDEFINED:
1277 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(),
deferred_logger );
1296 const auto current =
ws.production_cmode;
1297 const auto&
controls = well.productionControls(summaryState);
1299 case Well::ProducerCMode::ORAT:
1305 for (
int p = 0;
p<
np; ++
p) {
1309 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1312 for (
int p = 0;
p<
np; ++
p) {
1319 case Well::ProducerCMode::WRAT:
1325 for (
int p = 0;
p<
np; ++
p) {
1329 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1332 for (
int p = 0;
p<
np; ++
p) {
1339 case Well::ProducerCMode::GRAT:
1345 for (
int p = 0;
p<
np; ++
p) {
1349 const std::vector<Scalar >
fractions = initialWellRateFractions(simulator, well_state);
1352 for (
int p = 0;
p<
np; ++
p) {
1361 case Well::ProducerCMode::LRAT:
1364 -
ws.surface_rates[ pu.phase_pos[Oil] ];
1368 for (
int p = 0;
p<
np; ++
p) {
1372 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1375 for (
int p = 0;
p<
np; ++
p) {
1382 case Well::ProducerCMode::CRAT:
1384 OPM_DEFLOG_THROW(std::runtime_error,
1385 fmt::format(
"CRAT control not supported, well {}", this->name()),
1388 case Well::ProducerCMode::RESV:
1390 std::vector<Scalar>
convert_coeff(this->number_of_phases_, 1.0);
1391 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_,
ws.surface_rates,
convert_coeff);
1393 for (
int p = 0;
p<
np; ++
p) {
1400 for (
int p = 0;
p<
np; ++
p) {
1404 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1405 for (
int p = 0;
p<
np; ++
p) {
1410 std::vector<Scalar>
hrates(this->number_of_phases_,0.);
1411 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1414 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1417 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1420 std::vector<Scalar>
hrates_resv(this->number_of_phases_,0.);
1421 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_,
hrates,
hrates_resv);
1426 for (
int p = 0;
p<
np; ++
p) {
1430 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1431 for (
int p = 0;
p<
np; ++
p) {
1438 case Well::ProducerCMode::BHP:
1442 for (
int p = 0;
p<
np; ++
p) {
1449 for (
int p = 0;
p<
np; ++
p) {
1450 ws.surface_rates[
p] = -
ws.well_potentials[
p];
1455 case Well::ProducerCMode::THP:
1463 auto rates =
ws.surface_rates;
1464 this->adaptRatesForVFP(rates);
1465 const Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1468 ws.thp = this->getTHPConstraint(summaryState);
1472 const Scalar
total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1474 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1475 ws.surface_rates[
p] = -
ws.well_potentials[
p];
1481 case Well::ProducerCMode::GRUP:
1483 assert(well.isAvailableForGroupControl());
1484 const auto& group = schedule.getGroup(well.groupName(),
this->currentStep());
1486 well_state[well.name()].efficiency_scaling_factor;
1487 Scalar scale = this->getGroupProductionTargetRate(group,
1497 for (
int p = 0;
p<
np; ++
p) {
1498 ws.surface_rates[
p] *= scale;
1500 ws.trivial_group_target =
false;
1504 ws.trivial_group_target =
true;
1508 case Well::ProducerCMode::CMODE_UNDEFINED:
1509 case Well::ProducerCMode::NONE:
1511 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() ,
deferred_logger);
1522 template<
typename TypeTag>
1524 WellInterface<TypeTag>::
1525 wellUnderZeroRateTarget(
const Simulator& simulator,
1534 const auto& summaryState = simulator.vanguard().summaryState();
1535 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1541 template <
typename TypeTag>
1543 WellInterface<TypeTag>::wellUnderZeroGroupRateTarget(
const Simulator& simulator,
1551 const auto& summaryState = simulator.vanguard().summaryState();
1552 const auto& group_state = simulator.problem().wellModel().groupState();
1553 const auto& schedule = simulator.vanguard().schedule();
1554 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state,
deferred_logger);
1559 template<
typename TypeTag>
1561 WellInterface<TypeTag>::
1562 stoppedOrZeroRateTarget(
const Simulator& simulator,
1568 return this->wellIsStopped()
1569 || this->wellUnderZeroRateTarget(simulator, well_state,
deferred_logger);
1572 template<
typename TypeTag>
1573 std::vector<typename WellInterface<TypeTag>::Scalar>
1574 WellInterface<TypeTag>::
1575 initialWellRateFractions(
const Simulator& simulator,
1579 const int np = this->number_of_phases_;
1581 const auto&
ws = well_state.well(this->index_of_well_);
1584 for (
int p = 0;
p<
np; ++
p) {
1588 for (
int p = 0;
p<
np; ++
p) {
1596 const int nperf = this->number_of_local_perforations_;
1604 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1605 const auto& fs = intQuants.fluidState();
1608 for (
int p = 0;
p <
np; ++
p) {
1612 for (
int p = 0;
p <
np; ++
p) {
1622 template <
typename TypeTag>
1632 auto&
ws = well_state.well(this->index_of_well_);
1635 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1652 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1653 ws.surface_rates[
p] =
well_q_s[this->flowPhaseToModelCompIdx(
p)];
1662 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1664 const int comp_idx = this->flowPhaseToModelCompIdx(
p);
1665 Scalar& rate =
ws.surface_rates[
p];
1672 template <
typename TypeTag>
1673 std::vector<typename WellInterface<TypeTag>::Scalar>
1676 const IntensiveQuantities& intQuants,
1683 if (
static_cast<std::size_t
>(
perf) >= this->well_cells_.size()) {
1684 OPM_THROW(std::invalid_argument,
"The perforation index exceeds the size of the local containers - possibly wellIndex was called with a global instead of a local perforation index!");
1686 auto wi = std::vector<Scalar>
1689 if constexpr (! Indices::gasEnabled) {
1693 const auto&
wdfac = this->well_ecl_.getWDFAC();
1695 if (!
wdfac.useDFactor() || (this->well_index_[
perf] == 0.0)) {
1699 const Scalar
d = this->computeConnectionDFactor(
perf, intQuants,
ws);
1706 const auto&
connection = this->well_ecl_.getConnections()[
ws.perf_data.ecl_index[
perf]];
1709 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1714 const Scalar
invB =
getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1715 const Scalar
mob_g =
getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) *
invB;
1722 const Scalar
r2n =
b*
b + 4*
a*
c;
1724 const Scalar
rn = std::sqrt(
r2n);
1725 const Scalar
xn1 = (
b-
rn)*0.5/
a;
1729 const Scalar
xn2 = (
b+
rn)*0.5/
a;
1736 const Scalar
r2p =
b*
b - 4*
a*
c;
1738 const Scalar
rp = std::sqrt(
r2p);
1739 const Scalar
xp1 = (
rp-
b)*0.5/
a;
1743 const Scalar
xp2 = -(
rp+
b)*0.5/
a;
1753 template <
typename TypeTag>
1755 WellInterface<TypeTag>::
1756 updateConnectionDFactor(
const Simulator& simulator,
1759 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1763 auto&
d_factor =
ws.perf_data.connection_d_factor;
1765 for (
int perf = 0;
perf < this->number_of_local_perforations_; ++
perf) {
1767 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1773 template <
typename TypeTag>
1774 typename WellInterface<TypeTag>::Scalar
1775 WellInterface<TypeTag>::
1776 computeConnectionDFactor(
const int perf,
1777 const IntensiveQuantities& intQuants,
1781 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx,
regIdx);
1786 temperature =
ws.temperature,
1787 regIdx = this->pvtRegionIdx(), &intQuants]()
1789 const auto rv =
getValue(intQuants.fluidState().Rv());
1791 const auto&
gasPvt = FluidSystem::gasPvt();
1796 const Scalar
rv_sat =
gasPvt.saturatedOilVaporizationFactor
1805 rv,
getValue(intQuants.fluidState().Rvw()));
1808 const auto&
connection = this->well_ecl_.getConnections()
1809 [
ws.perf_data.ecl_index[
perf]];
1815 template <
typename TypeTag>
1817 WellInterface<TypeTag>::
1818 updateConnectionTransmissibilityFactor(
const Simulator& simulator,
1822 &
conns = this->well_ecl_.getConnections()]
1828 auto&
tmult =
ws.perf_data.connection_compaction_tmult;
1829 auto& ctf =
ws.perf_data.connection_transmissibility_factor;
1831 for (
int perf = 0;
perf < this->number_of_local_perforations_; ++
perf) {
1834 const auto& intQuants = simulator.model()
1845 template<
typename TypeTag>
1846 typename WellInterface<TypeTag>::Eval
1847 WellInterface<TypeTag>::getPerfCellPressure(
const typename WellInterface<TypeTag>::FluidState& fs)
const
1849 if constexpr (Indices::oilEnabled) {
1850 return fs.pressure(FluidSystem::oilPhaseIdx);
1851 }
else if constexpr (Indices::gasEnabled) {
1852 return fs.pressure(FluidSystem::gasPhaseIdx);
1854 return fs.pressure(FluidSystem::waterPhaseIdx);
1858 template <
typename TypeTag>
1859 template<
class Value,
class Callback>
1861 WellInterface<TypeTag>::
1862 getMobility(
const Simulator& simulator,
1864 std::vector<Value>&
mob,
1870 if constexpr (std::is_same_v<Value, Scalar>) {
1871 return std::array<Scalar,3>{};
1873 return std::array<Eval,3>{};
1876 if (
static_cast<std::size_t
>(
perf) >= this->well_cells_.size()) {
1877 OPM_THROW(std::invalid_argument,
"The perforation index exceeds the size of the local containers - possibly getMobility was called with a global instead of a local perforation index!");
1881 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1882 const auto& materialLawManager = simulator.problem().materialLawManager();
1886 const int satid = this->saturation_table_number_[
perf] - 1;
1890 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
1894 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
1897 if constexpr (has_solvent) {
1898 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1910 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
1914 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
1919 if constexpr (has_solvent) {
1920 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent",
deferred_logger);
1924 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1926 const auto&
connections = this->well_ecl_.getConnections();
1930 val *= this->inj_fc_multiplier_[
perf];
1937 template<
typename TypeTag>
1939 WellInterface<TypeTag>::
1940 updateWellStateWithTHPTargetProd(
const Simulator& simulator,
1945 const auto&
summary_state = simulator.vanguard().summaryState();
1950 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
1951 if (thp_update_iterations) {
1958 auto&
ws = well_state.well(this->name());
1959 ws.surface_rates = rates;
1968 template <
typename TypeTag>
1970 WellInterface<TypeTag>::
1971 computeConnLevelProdInd(
const FluidState& fs,
1972 const std::function<Scalar(
const Scalar)>&
connPICalc,
1973 const std::vector<Scalar>& mobility,
1977 const int np = this->number_of_phases_;
1978 for (
int p = 0;
p <
np; ++
p) {
1982 mobility[this->flowPhaseToModelCompIdx(
p)]
1983 * fs.invB(this->flowPhaseToModelPhaseIdx(
p)).value();
1988 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1989 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1991 const auto io = pu.phase_pos[Oil];
1992 const auto ig = pu.phase_pos[Gas];
2003 template <
typename TypeTag>
2005 WellInterface<TypeTag>::
2006 computeConnLevelInjInd(
const FluidState& fs,
2008 const std::function<Scalar(
const Scalar)>&
connIICalc,
2009 const std::vector<Scalar>& mobility,
2018 phase_pos = pu.phase_pos[Gas];
2021 phase_pos = pu.phase_pos[Oil];
2024 phase_pos = pu.phase_pos[Water];
2028 fmt::format(
"Unsupported Injector Type ({}) "
2029 "for well {} during connection I.I. calculation",
2034 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2038 template<
typename TypeTag>
2039 template<
class GasLiftSingleWell>
2040 std::unique_ptr<GasLiftSingleWell>
2041 WellInterface<TypeTag>::
2042 initializeGliftWellTest_(
const Simulator& simulator,
2050 auto& comm = simulator.vanguard().grid().comm();
2051 ecl_well_map.try_emplace(this->name(), &(this->wellEcl()), this->indexOfWell());
2054 simulator.vanguard().schedule(),
2055 simulator.vanguard().summaryState(),
2056 simulator.episodeIndex(),
2057 simulator.model().newtonMethod().numIterations(),
2068 const auto&
summary_state = simulator.vanguard().summaryState();
2069 return std::make_unique<GasLiftSingleWell>(*
this,
Definition DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition ParallelWellInfo.hpp:195
Definition SingleWellState.hpp:42
Definition WellInterfaceIndices.hpp:34
Definition WellInterface.hpp:77
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:58
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:1625
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:66
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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Static data associated with a well perforation.
Definition PerforationData.hpp:30
Definition BlackoilPhases.hpp:46