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_();
176 template<
typename TypeTag>
177 typename WellInterface<TypeTag>::Scalar
178 WellInterface<TypeTag>::
181 if constexpr (has_micp) {
182 return this->wurea_();
188 template<
typename TypeTag>
190 WellInterface<TypeTag>::
191 updateWellControl(
const Simulator& simulator,
192 const IndividualOrGroup
iog,
202 const auto& summaryState = simulator.vanguard().summaryState();
203 const auto& schedule = simulator.vanguard().schedule();
204 const auto& well = this->well_ecl_;
205 auto&
ws = well_state.well(this->index_of_well_);
207 if (well.isInjector()) {
212 bool oscillating = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) >= this->param_.max_number_of_well_switches_;
213 const int episodeIdx = simulator.episodeIndex();
214 const int iterationIdx = simulator.model().newtonMethod().numIterations();
218 bool output = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) == this->param_.max_number_of_well_switches_;
220 std::ostringstream
ss;
221 ss <<
" The control mode for well " << this->name()
222 <<
" is oscillating\n"
223 <<
" We don't allow for more than "
224 << this->param_.max_number_of_well_switches_
225 <<
" switches. The control is kept at " <<
from;
228 this->well_control_log_.push_back(
from);
233 if (
iog == IndividualOrGroup::Individual) {
235 }
else if (
iog == IndividualOrGroup::Group) {
241 Parallel::Communication
cc = simulator.vanguard().grid().comm();
245 if (well.isInjector()) {
250 std::ostringstream
ss;
251 ss <<
" Switching control mode for well " << this->name()
255 ss <<
" on rank " <<
cc.rank();
259 this->well_control_log_.push_back(
from);
260 updateWellStateWithTarget(simulator, group_state, well_state,
deferred_logger);
267 template<
typename TypeTag>
269 WellInterface<TypeTag>::
270 updateWellControlAndStatusLocalIteration(
const Simulator& simulator,
281 const auto&
summary_state = simulator.vanguard().summaryState();
282 const auto& schedule = simulator.vanguard().schedule();
283 auto&
ws = well_state.well(this->index_of_well_);
285 if (this->isInjector()) {
290 const bool oscillating = std::count(this->well_control_log_.begin(),
this->well_control_log_.end(),
from) >= this->param_.max_number_of_well_switches_;
292 if (
oscillating || this->wellUnderZeroRateTarget(simulator, well_state,
deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
296 const Scalar
sgn = this->isInjector() ? 1.0 : -1.0;
297 if (!this->wellIsStopped()){
310 bool isGroupControl =
ws.production_cmode == Well::ProducerCMode::GRUP ||
ws.injection_cmode == Well::InjectorCMode::GRUP;
311 if (! (
isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
314 if (
hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
319 const bool thp_controlled = this->isInjector() ?
ws.injection_cmode == Well::InjectorCMode::THP :
320 ws.production_cmode == Well::ProducerCMode::THP;
323 updateWellStateWithTarget(simulator, group_state, well_state,
deferred_logger);
334 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
339 std::vector<Scalar> rates(this->num_components_);
340 if (this->isInjector()){
341 const Scalar
bhp_thp = WellBhpThpCalculator(*this).
342 calculateBhpFromThp(well_state, rates,
345 this->getRefDensity(),
351 const Scalar
bhp_min = WellBhpThpCalculator(*this).
352 calculateMinimumBhpFromThp(well_state,
355 this->getRefDensity());
364 well_state.well(this->index_of_well_).thp = this->getTHPConstraint(
summary_state);
375 template<
typename TypeTag>
377 WellInterface<TypeTag>::
378 wellTesting(
const Simulator& simulator,
394 const auto&
summary_state = simulator.vanguard().summaryState();
397 ws.production_cmode = Well::ProducerCMode::THP;
400 ws.production_cmode = Well::ProducerCMode::BHP;
407 if (this->isProducer()) {
408 const auto& schedule = simulator.vanguard().schedule();
409 const auto report_step = simulator.episodeIndex();
410 const auto&
glo = schedule.glo(report_step);
412 gliftBeginTimeStepWellTestUpdateALQ(simulator,
431 const auto msg = fmt::format(
"WTEST: Well {} is not solvable (physical)", this->name());
438 if ( !this->isOperableAndSolvable() ) {
439 const auto msg = fmt::format(
"WTEST: Well {} is not operable (physical)", this->name());
446 }
catch (
const std::exception&
e) {
447 const std::string
msg = fmt::format(
"well {}: computeWellPotentials() "
448 "failed during testing for re-opening: ",
449 this->name(),
e.what());
454 for (
int p = 0;
p <
np; ++
p) {
478 well_test_state.open_well(this->name());
480 std::string
msg = std::string(
"well ") + this->name() + std::string(
" is re-opened");
486 well_test_state.open_completion(this->name(),
completion.first);
491 open_times.try_emplace(this->name(), well_test_state.lastTestTime(
this->name()));
498 template<
typename TypeTag>
500 WellInterface<TypeTag>::
501 iterateWellEquations(
const Simulator& simulator,
508 const auto&
summary_state = simulator.vanguard().summaryState();
509 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(
summary_state) : Well::InjectionControls(0);
510 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(
summary_state) : Well::ProductionControls(0);
511 bool converged =
false;
514 if (!this->param_.local_well_solver_control_switching_){
517 if (this->param_.use_implicit_ipr_ &&
this->well_ecl_.isProducer() &&
this->wellHasTHPConstraints(
summary_state) && (
this->well_ecl_.getStatus() == WellStatus::OPEN)) {
525 const std::string
msg =
"Inner well iterations failed for well " + this->name() +
" Treat the well as unconverged. ";
532 template<
typename TypeTag>
534 WellInterface<TypeTag>::
535 solveWellWithTHPConstraint(
const Simulator& simulator,
544 const auto&
summary_state = simulator.vanguard().summaryState();
546 bool converged =
true;
547 auto&
ws = well_state.well(this->index_of_well_);
549 if (this->wellIsStopped()) {
554 const auto msg = fmt::format(
"estimateOperableBhp: Did not find operable BHP for well {}", this->name());
563 const Scalar bhp = std::max(
bhp_target.value(),
573 const bool isThp =
ws.production_cmode == Well::ProducerCMode::THP;
576 auto rates = well_state.well(this->index_of_well_).surface_rates;
577 this->adaptRatesForVFP(rates);
582 this->operability_status_.use_vfpexplicit =
true;
585 const Scalar reltol = 1
e-3;
588 const auto msg = fmt::format(
"Well {} converged to an unstable solution, re-solving", this->name());
600 this->operability_status_.use_vfpexplicit =
true;
611 const Scalar bhp = std::max(
bhp_target.value(),
615 converged = this->iterateWellEqWithSwitching(simulator,
dt,
625 this->operability_status_.can_obtain_bhp_with_thp_limit =
is_operable;
626 this->operability_status_.obey_thp_limit_under_bhp_limit =
is_operable;
630 template<
typename TypeTag>
631 std::optional<typename WellInterface<TypeTag>::Scalar>
632 WellInterface<TypeTag>::
633 estimateOperableBhp(
const Simulator& simulator,
642 Scalar
bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state,
this->well_ecl_,
summary_state,
this->getRefDensity());
645 if (!converged || this->wellIsStopped()) {
649 auto rates = well_state.well(this->index_of_well_).surface_rates;
650 this->adaptRatesForVFP(rates);
651 return WellBhpThpCalculator(*this).estimateStableBhp(well_state,
this->well_ecl_, rates,
this->getRefDensity(),
summary_state);
654 template<
typename TypeTag>
656 WellInterface<TypeTag>::
657 solveWellWithBhp(
const Simulator& simulator,
668 auto&
ws = well_state.well(this->index_of_well_);
671 if (this->isInjector()) {
675 ws.injection_cmode = Well::InjectorCMode::BHP;
680 ws.production_cmode = Well::ProducerCMode::BHP;
691 template<
typename TypeTag>
693 WellInterface<TypeTag>::
694 solveWellWithZeroRate(
const Simulator& simulator,
712 template<
typename TypeTag>
714 WellInterface<TypeTag>::
715 solveWellForTesting(
const Simulator& simulator,
721 const double dt = simulator.timeStepSize();
722 const bool converged = iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
724 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" converged");
727 const int max_iter = this->param_.max_welleq_iter_;
728 deferred_logger.debug(
"WellTest: Well equation for well " + this->name() +
" failed converging in "
729 + std::to_string(
max_iter) +
" iterations");
734 template<
typename TypeTag>
736 WellInterface<TypeTag>::
737 solveWellEquation(
const Simulator& simulator,
743 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
748 const double dt = simulator.timeStepSize();
749 bool converged = iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
759 auto&
ws = well_state.well(this->indexOfWell());
761 if (this->well_ecl_.isInjector()) {
762 thp_control =
ws.injection_cmode == Well::InjectorCMode::THP;
764 ws.injection_cmode = Well::InjectorCMode::BHP;
768 thp_control =
ws.production_cmode == Well::ProducerCMode::THP;
770 ws.production_cmode = Well::ProducerCMode::BHP;
775 const std::string
msg = std::string(
"The newly opened well ") + this->name()
776 + std::string(
" with THP control did not converge during inner iterations, we try again with bhp control");
778 converged = this->iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
783 const int max_iter = this->param_.max_welleq_iter_;
784 deferred_logger.debug(
"Compute initial well solution for well " + this->name() +
". Failed to converge in "
785 + std::to_string(
max_iter) +
" iterations");
792 template <
typename TypeTag>
794 WellInterface<TypeTag>::
795 assembleWellEq(
const Simulator& simulator,
802 prepareWellBeforeAssembling(simulator,
dt, well_state, group_state,
deferred_logger);
803 assembleWellEqWithoutIteration(simulator,
dt, well_state, group_state,
deferred_logger);
808 template <
typename TypeTag>
810 WellInterface<TypeTag>::
811 assembleWellEqWithoutIteration(
const Simulator& simulator,
818 const auto&
summary_state = simulator.vanguard().summaryState();
819 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(
summary_state) : Well::InjectionControls(0);
820 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(
summary_state) : Well::ProductionControls(0);
828 template<
typename TypeTag>
830 WellInterface<TypeTag>::
831 prepareWellBeforeAssembling(
const Simulator& simulator,
838 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
840 if (this->param_.check_well_operability_iter_)
844 const int iteration_idx = simulator.model().newtonMethod().numIterations();
846 const auto&
ws = well_state.well(this->indexOfWell());
851 std::any_of(
ws.surface_rates.begin(),
852 ws.surface_rates.begin() + well_state.numPhases(),
853 [](Scalar rate) { return rate != Scalar(0.0); });
855 this->operability_status_.solvable =
true;
856 bool converged = this->iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
864 this->operability_status_.resetOperability();
866 deferred_logger.debug(
" " + this->name() +
" is re-opened after being stopped during local solve");
871 if (this->isInjector()) {
878 deferred_logger.debug(
" " + this->name() +
" switched from " +
from +
" to " +
to +
" during local solve");
883 if (this->param_.shut_unsolvable_wells_)
884 this->operability_status_.solvable =
false;
887 if (this->operability_status_.has_negative_potentials) {
892 }
catch (
const std::exception&
e) {
893 const std::string
msg = fmt::format(
"well {}: computeWellPotentials() failed "
894 "during attempt to recompute potentials for well: ",
895 this->name(),
e.what());
897 this->operability_status_.has_negative_potentials =
true;
899 auto&
ws = well_state.well(this->indexOfWell());
900 const int np = well_state.numPhases();
901 for (
int p = 0;
p <
np; ++
p) {
905 this->changed_to_open_this_step_ =
false;
906 const bool well_operable = this->operability_status_.isOperableAndSolvable();
909 deferred_logger.debug(
" well " + this->name() +
" gets STOPPED during iteration ");
911 changed_to_stopped_this_step_ =
true;
913 deferred_logger.debug(
" well " + this->name() +
" gets REVIVED during iteration ");
915 changed_to_stopped_this_step_ =
false;
916 this->changed_to_open_this_step_ =
true;
920 template<
typename TypeTag>
922 WellInterface<TypeTag>::addCellRates(RateVector& rates,
int cellIdx)
const
924 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
929 for (
int i = 0; i < RateVector::dimension; ++i) {
930 rates[i] += connectionRates_[
perfIdx][i];
936 template<
typename TypeTag>
937 typename WellInterface<TypeTag>::Scalar
938 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(
int cellIdx,
int phaseIdx)
const
942 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
947 OPM_THROW(std::invalid_argument,
"The well with name " + this->name()
948 +
" does not perforate cell " + std::to_string(
cellIdx));
955 template<
typename TypeTag>
957 WellInterface<TypeTag>::
958 checkWellOperability(
const Simulator& simulator,
963 if (!this->param_.check_well_operability_) {
967 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
972 if (!this->operability_status_.isOperableAndSolvable()) {
973 this->operability_status_.use_vfpexplicit =
true;
975 "well not operable, trying with explicit vfp lookup: " + this->name());
980 template<
typename TypeTag>
982 WellInterface<TypeTag>::
983 gliftBeginTimeStepWellTestIterateWellEquations(
const Simulator& simulator,
990 const auto& well_name = this->name();
991 assert(this->wellHasTHPConstraints(simulator.vanguard().summaryState()));
992 const auto& schedule = simulator.vanguard().schedule();
995 if(
glo.active() &&
glo.has_well(well_name)) {
996 const auto increment =
glo.gaslift_increment();
997 auto alq = well_state.getALQ(well_name);
1000 well_state.setALQ(well_name, alq);
1011 return iterateWellEquations(simulator,
dt, well_state, group_state,
deferred_logger);
1015 template<
typename TypeTag>
1017 WellInterface<TypeTag>::
1018 gliftBeginTimeStepWellTestUpdateALQ(
const Simulator& simulator,
1026 const auto&
summary_state = simulator.vanguard().summaryState();
1027 const auto& well_name = this->name();
1029 const std::string
msg = fmt::format(
"GLIFT WTEST: Well {} does not have THP constraints", well_name);
1033 const auto& schedule = simulator.vanguard().schedule();
1036 if (!
glo.has_well(well_name)) {
1037 const std::string
msg = fmt::format(
1038 "GLIFT WTEST: Well {} : Gas lift not activated: "
1039 "WLIFTOPT is probably missing. Skipping.", well_name);
1046 std::unique_ptr<GasLiftSingleWell>
glift =
1057 well_state.setALQ(well_name,
wtest_alq);
1059 "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1065 "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1067 unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.getALQ(well_name)));
1072 "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1079 template<
typename TypeTag>
1081 WellInterface<TypeTag>::
1082 updateWellOperability(
const Simulator& simulator,
1087 if (this->param_.local_well_solver_control_switching_) {
1088 const bool success = updateWellOperabilityFromWellEq(simulator, well_state,
deferred_logger);
1092 deferred_logger.debug(
"Operability check using well equations did not converge for well "
1093 + this->name() +
", reverting to classical approach." );
1096 this->operability_status_.resetOperability();
1098 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1099 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1100 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1101 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1108 checkOperabilityUnderBHPLimit(well_state, simulator,
deferred_logger);
1112 checkOperabilityUnderTHPLimit(simulator, well_state,
deferred_logger);
1116 template<
typename TypeTag>
1118 WellInterface<TypeTag>::
1119 updateWellOperabilityFromWellEq(
const Simulator& simulator,
1125 assert(this->param_.local_well_solver_control_switching_);
1126 this->operability_status_.resetOperability();
1128 const auto& group_state = simulator.problem().wellModel().groupState();
1129 const double dt = simulator.timeStepSize();
1135 template<
typename TypeTag>
1137 WellInterface<TypeTag>::
1138 updateWellStateWithTarget(
const Simulator& simulator,
1145 const auto& well = this->well_ecl_;
1146 const int well_index = this->index_of_well_;
1147 auto&
ws = well_state.well(well_index);
1149 const int np = well_state.numPhases();
1150 const auto& summaryState = simulator.vanguard().summaryState();
1151 const auto& schedule = simulator.vanguard().schedule();
1153 if (this->wellIsStopped()) {
1154 for (
int p = 0;
p<
np; ++
p) {
1155 ws.surface_rates[
p] = 0;
1161 if (this->isInjector() )
1163 const auto&
controls = well.injectionControls(summaryState);
1168 case InjectorType::WATER:
1170 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
1173 case InjectorType::OIL:
1175 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
1178 case InjectorType::GAS:
1180 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
1184 OPM_DEFLOG_THROW(std::runtime_error,
"Expected WATER, OIL or GAS as type for injectors " + this->name(),
deferred_logger );
1187 const auto current =
ws.injection_cmode;
1190 case Well::InjectorCMode::RATE:
1193 if(this->rsRvInj() > 0) {
1194 if (
injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1195 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] =
controls.surface_rate * this->rsRvInj();
1196 }
else if (
injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1197 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] =
controls.surface_rate * this->rsRvInj();
1199 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 );
1205 case Well::InjectorCMode::RESV:
1207 std::vector<Scalar>
convert_coeff(this->number_of_phases_, 1.0);
1208 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_,
convert_coeff);
1214 case Well::InjectorCMode::THP:
1216 auto rates =
ws.surface_rates;
1217 Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
1221 this->getRefDensity(),
1224 ws.thp = this->getTHPConstraint(summaryState);
1229 Scalar
total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1231 ws.surface_rates =
ws.well_potentials;
1235 case Well::InjectorCMode::BHP:
1239 for (
int p = 0;
p<
np; ++
p) {
1246 ws.surface_rates =
ws.well_potentials;
1250 case Well::InjectorCMode::GRUP:
1252 assert(well.isAvailableForGroupControl());
1253 const auto& group = schedule.getGroup(well.groupName(),
this->currentStep());
1255 well_state[well.name()].efficiency_scaling_factor;
1256 std::optional<Scalar>
target =
1257 this->getGroupInjectionTargetRate(group,
1269 case Well::InjectorCMode::CMODE_UNDEFINED:
1271 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name(),
deferred_logger );
1290 const auto current =
ws.production_cmode;
1291 const auto&
controls = well.productionControls(summaryState);
1293 case Well::ProducerCMode::ORAT:
1299 for (
int p = 0;
p<
np; ++
p) {
1303 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1306 for (
int p = 0;
p<
np; ++
p) {
1313 case Well::ProducerCMode::WRAT:
1319 for (
int p = 0;
p<
np; ++
p) {
1323 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1326 for (
int p = 0;
p<
np; ++
p) {
1333 case Well::ProducerCMode::GRAT:
1339 for (
int p = 0;
p<
np; ++
p) {
1343 const std::vector<Scalar >
fractions = initialWellRateFractions(simulator, well_state);
1346 for (
int p = 0;
p<
np; ++
p) {
1355 case Well::ProducerCMode::LRAT:
1358 -
ws.surface_rates[ pu.phase_pos[Oil] ];
1362 for (
int p = 0;
p<
np; ++
p) {
1366 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1369 for (
int p = 0;
p<
np; ++
p) {
1376 case Well::ProducerCMode::CRAT:
1378 OPM_DEFLOG_THROW(std::runtime_error,
1379 fmt::format(
"CRAT control not supported, well {}", this->name()),
1382 case Well::ProducerCMode::RESV:
1384 std::vector<Scalar>
convert_coeff(this->number_of_phases_, 1.0);
1385 this->rateConverter_.calcCoeff( 0, this->pvtRegionIdx_,
ws.surface_rates,
convert_coeff);
1387 for (
int p = 0;
p<
np; ++
p) {
1394 for (
int p = 0;
p<
np; ++
p) {
1398 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1399 for (
int p = 0;
p<
np; ++
p) {
1404 std::vector<Scalar>
hrates(this->number_of_phases_,0.);
1405 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1408 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1411 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1414 std::vector<Scalar>
hrates_resv(this->number_of_phases_,0.);
1415 this->rateConverter_.calcReservoirVoidageRates( 0, this->pvtRegionIdx_,
hrates,
hrates_resv);
1420 for (
int p = 0;
p<
np; ++
p) {
1424 const std::vector<Scalar>
fractions = initialWellRateFractions(simulator, well_state);
1425 for (
int p = 0;
p<
np; ++
p) {
1432 case Well::ProducerCMode::BHP:
1436 for (
int p = 0;
p<
np; ++
p) {
1443 for (
int p = 0;
p<
np; ++
p) {
1444 ws.surface_rates[
p] = -
ws.well_potentials[
p];
1449 case Well::ProducerCMode::THP:
1457 auto rates =
ws.surface_rates;
1458 this->adaptRatesForVFP(rates);
1459 const Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1462 ws.thp = this->getTHPConstraint(summaryState);
1466 const Scalar
total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1468 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1469 ws.surface_rates[
p] = -
ws.well_potentials[
p];
1475 case Well::ProducerCMode::GRUP:
1477 assert(well.isAvailableForGroupControl());
1478 const auto& group = schedule.getGroup(well.groupName(),
this->currentStep());
1480 well_state[well.name()].efficiency_scaling_factor;
1481 Scalar scale = this->getGroupProductionTargetRate(group,
1491 for (
int p = 0;
p<
np; ++
p) {
1492 ws.surface_rates[
p] *= scale;
1494 ws.trivial_group_target =
false;
1498 ws.trivial_group_target =
true;
1502 case Well::ProducerCMode::CMODE_UNDEFINED:
1503 case Well::ProducerCMode::NONE:
1505 OPM_DEFLOG_THROW(std::runtime_error,
"Well control must be specified for well " + this->name() ,
deferred_logger);
1516 template<
typename TypeTag>
1518 WellInterface<TypeTag>::
1519 wellUnderZeroRateTarget(
const Simulator& simulator,
1528 const auto& summaryState = simulator.vanguard().summaryState();
1529 return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1535 template <
typename TypeTag>
1537 WellInterface<TypeTag>::wellUnderZeroGroupRateTarget(
const Simulator& simulator,
1545 const auto& summaryState = simulator.vanguard().summaryState();
1546 const auto& group_state = simulator.problem().wellModel().groupState();
1547 const auto& schedule = simulator.vanguard().schedule();
1548 return this->zeroGroupRateTarget(summaryState, schedule, well_state, group_state,
deferred_logger);
1553 template<
typename TypeTag>
1555 WellInterface<TypeTag>::
1556 stoppedOrZeroRateTarget(
const Simulator& simulator,
1562 return this->wellIsStopped()
1563 || this->wellUnderZeroRateTarget(simulator, well_state,
deferred_logger);
1566 template<
typename TypeTag>
1567 std::vector<typename WellInterface<TypeTag>::Scalar>
1568 WellInterface<TypeTag>::
1569 initialWellRateFractions(
const Simulator& simulator,
1573 const int np = this->number_of_phases_;
1575 const auto&
ws = well_state.well(this->index_of_well_);
1578 for (
int p = 0;
p<
np; ++
p) {
1582 for (
int p = 0;
p<
np; ++
p) {
1590 const int nperf = this->number_of_local_perforations_;
1598 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1599 const auto& fs = intQuants.fluidState();
1602 for (
int p = 0;
p <
np; ++
p) {
1606 for (
int p = 0;
p <
np; ++
p) {
1616 template <
typename TypeTag>
1626 auto&
ws = well_state.well(this->index_of_well_);
1629 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1646 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1647 ws.surface_rates[
p] =
well_q_s[this->flowPhaseToModelCompIdx(
p)];
1656 for (
int p = 0;
p < this->number_of_phases_; ++
p) {
1658 const int comp_idx = this->flowPhaseToModelCompIdx(
p);
1659 Scalar& rate =
ws.surface_rates[
p];
1666 template <
typename TypeTag>
1667 std::vector<typename WellInterface<TypeTag>::Scalar>
1670 const IntensiveQuantities& intQuants,
1677 if (
static_cast<std::size_t
>(
perf) >= this->well_cells_.size()) {
1678 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!");
1680 auto wi = std::vector<Scalar>
1683 if constexpr (! Indices::gasEnabled) {
1687 const auto&
wdfac = this->well_ecl_.getWDFAC();
1689 if (!
wdfac.useDFactor() || (this->well_index_[
perf] == 0.0)) {
1693 const Scalar
d = this->computeConnectionDFactor(
perf, intQuants,
ws);
1700 const auto&
connection = this->well_ecl_.getConnections()[
ws.perf_data.ecl_index[
perf]];
1703 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1708 const Scalar
invB =
getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1709 const Scalar
mob_g =
getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) *
invB;
1716 const Scalar
r2n =
b*
b + 4*
a*
c;
1718 const Scalar
rn = std::sqrt(
r2n);
1719 const Scalar
xn1 = (
b-
rn)*0.5/
a;
1723 const Scalar
xn2 = (
b+
rn)*0.5/
a;
1730 const Scalar
r2p =
b*
b - 4*
a*
c;
1732 const Scalar
rp = std::sqrt(
r2p);
1733 const Scalar
xp1 = (
rp-
b)*0.5/
a;
1737 const Scalar
xp2 = -(
rp+
b)*0.5/
a;
1747 template <
typename TypeTag>
1749 WellInterface<TypeTag>::
1750 updateConnectionDFactor(
const Simulator& simulator,
1753 if (! this->well_ecl_.getWDFAC().useDFactor()) {
1757 auto&
d_factor =
ws.perf_data.connection_d_factor;
1759 for (
int perf = 0;
perf < this->number_of_local_perforations_; ++
perf) {
1761 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1767 template <
typename TypeTag>
1768 typename WellInterface<TypeTag>::Scalar
1769 WellInterface<TypeTag>::
1770 computeConnectionDFactor(
const int perf,
1771 const IntensiveQuantities& intQuants,
1775 return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx,
regIdx);
1780 temperature =
ws.temperature,
1781 regIdx = this->pvtRegionIdx(), &intQuants]()
1783 const auto rv =
getValue(intQuants.fluidState().Rv());
1785 const auto&
gasPvt = FluidSystem::gasPvt();
1790 const Scalar
rv_sat =
gasPvt.saturatedOilVaporizationFactor
1799 rv,
getValue(intQuants.fluidState().Rvw()));
1802 const auto&
connection = this->well_ecl_.getConnections()
1803 [
ws.perf_data.ecl_index[
perf]];
1809 template <
typename TypeTag>
1811 WellInterface<TypeTag>::
1812 updateConnectionTransmissibilityFactor(
const Simulator& simulator,
1816 &
conns = this->well_ecl_.getConnections()]
1822 auto&
tmult =
ws.perf_data.connection_compaction_tmult;
1823 auto& ctf =
ws.perf_data.connection_transmissibility_factor;
1825 for (
int perf = 0;
perf < this->number_of_local_perforations_; ++
perf) {
1828 const auto& intQuants = simulator.model()
1839 template<
typename TypeTag>
1840 typename WellInterface<TypeTag>::Eval
1841 WellInterface<TypeTag>::getPerfCellPressure(
const typename WellInterface<TypeTag>::FluidState& fs)
const
1843 if constexpr (Indices::oilEnabled) {
1844 return fs.pressure(FluidSystem::oilPhaseIdx);
1845 }
else if constexpr (Indices::gasEnabled) {
1846 return fs.pressure(FluidSystem::gasPhaseIdx);
1848 return fs.pressure(FluidSystem::waterPhaseIdx);
1852 template <
typename TypeTag>
1853 template<
class Value,
class Callback>
1855 WellInterface<TypeTag>::
1856 getMobility(
const Simulator& simulator,
1858 std::vector<Value>&
mob,
1864 if constexpr (std::is_same_v<Value, Scalar>) {
1865 return std::array<Scalar,3>{};
1867 return std::array<Eval,3>{};
1870 if (
static_cast<std::size_t
>(
perf) >= this->well_cells_.size()) {
1871 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!");
1875 const auto& intQuants = simulator.model().intensiveQuantities(
cell_idx, 0);
1876 const auto& materialLawManager = simulator.problem().materialLawManager();
1880 const int satid = this->saturation_table_number_[
perf] - 1;
1884 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
1888 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
1891 if constexpr (has_solvent) {
1892 mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1904 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
1908 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
1913 if constexpr (has_solvent) {
1914 OPM_DEFLOG_THROW(std::runtime_error,
"individual mobility for wells does not work in combination with solvent",
deferred_logger);
1918 if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
1920 const auto&
connections = this->well_ecl_.getConnections();
1924 val *= this->inj_fc_multiplier_[
perf];
1931 template<
typename TypeTag>
1933 WellInterface<TypeTag>::
1934 updateWellStateWithTHPTargetProd(
const Simulator& simulator,
1939 const auto&
summary_state = simulator.vanguard().summaryState();
1944 std::vector<Scalar> rates(this->number_of_phases_, 0.0);
1945 if (thp_update_iterations) {
1952 auto&
ws = well_state.well(this->name());
1953 ws.surface_rates = rates;
1962 template <
typename TypeTag>
1964 WellInterface<TypeTag>::
1965 computeConnLevelProdInd(
const FluidState& fs,
1966 const std::function<Scalar(
const Scalar)>&
connPICalc,
1967 const std::vector<Scalar>& mobility,
1971 const int np = this->number_of_phases_;
1972 for (
int p = 0;
p <
np; ++
p) {
1976 mobility[this->flowPhaseToModelCompIdx(
p)]
1977 * fs.invB(this->flowPhaseToModelPhaseIdx(
p)).value();
1982 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1983 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1985 const auto io = pu.phase_pos[Oil];
1986 const auto ig = pu.phase_pos[Gas];
1997 template <
typename TypeTag>
1999 WellInterface<TypeTag>::
2000 computeConnLevelInjInd(
const FluidState& fs,
2002 const std::function<Scalar(
const Scalar)>&
connIICalc,
2003 const std::vector<Scalar>& mobility,
2012 phase_pos = pu.phase_pos[Gas];
2015 phase_pos = pu.phase_pos[Oil];
2018 phase_pos = pu.phase_pos[Water];
2022 fmt::format(
"Unsupported Injector Type ({}) "
2023 "for well {} during connection I.I. calculation",
2028 const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2032 template<
typename TypeTag>
2033 template<
class GasLiftSingleWell>
2034 std::unique_ptr<GasLiftSingleWell>
2035 WellInterface<TypeTag>::
2036 initializeGliftWellTest_(
const Simulator& simulator,
2044 auto& comm = simulator.vanguard().grid().comm();
2045 ecl_well_map.try_emplace(this->name(), &(this->wellEcl()), this->indexOfWell());
2048 simulator.vanguard().schedule(),
2049 simulator.vanguard().summaryState(),
2050 simulator.episodeIndex(),
2051 simulator.model().newtonMethod().numIterations(),
2062 const auto&
summary_state = simulator.vanguard().summaryState();
2063 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:1619
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