26 template<
typename TypeTag>
30 , current_step_(time_step)
34 OPM_THROW(std::invalid_argument,
"Null pointer of Well is used to construct WellInterface");
38 OPM_THROW(std::invalid_argument,
"Negtive time step is used to construct WellInterface");
42 OPM_THROW(std::invalid_argument,
"Null pointer of Wells is used to construct WellInterface");
45 const std::string& well_name = well->name();
49 for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
50 if (well_name == std::string(wells->name[index_well])) {
57 assert(index_well != wells->number_of_wells);
59 index_of_well_ = index_well;
60 well_type_ = wells->type[index_well];
61 number_of_phases_ = wells->number_of_phases;
65 comp_frac_.resize(number_of_phases_);
66 const int index_begin = index_well * number_of_phases_;
67 std::copy(wells->comp_frac + index_begin,
68 wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
71 well_controls_ = wells->ctrls[index_well];
73 ref_depth_ = wells->depth_ref[index_well];
77 const int perf_index_begin = wells->well_connpos[index_well];
78 const int perf_index_end = wells->well_connpos[index_well + 1];
79 number_of_perforations_ = perf_index_end - perf_index_begin;
80 first_perf_ = perf_index_begin;
82 well_cells_.resize(number_of_perforations_);
83 std::copy(wells->well_cells + perf_index_begin,
84 wells->well_cells + perf_index_end,
85 well_cells_.begin() );
87 well_index_.resize(number_of_perforations_);
88 std::copy(wells->WI + perf_index_begin,
89 wells->WI + perf_index_end,
90 well_index_.begin() );
92 saturation_table_number_.resize(number_of_perforations_);
93 std::copy(wells->sat_table_id + perf_index_begin,
94 wells->sat_table_id + perf_index_end,
95 saturation_table_number_.begin() );
98 well_efficiency_factor_ = 1.0;
105 template<
typename TypeTag>
108 init(
const PhaseUsage* phase_usage_arg,
109 const std::vector<bool>* active_arg,
110 const std::vector<double>& ,
111 const double gravity_arg,
114 phase_usage_ = phase_usage_arg;
115 active_ = active_arg;
116 gravity_ = gravity_arg;
123 template<
typename TypeTag>
125 WellInterface<TypeTag>::
126 setVFPProperties(
const VFPProperties* vfp_properties_arg)
128 vfp_properties_ = vfp_properties_arg;
135 template<
typename TypeTag>
140 return well_ecl_->name();
147 template<
typename TypeTag>
159 template<
typename TypeTag>
164 return well_controls_;
171 template<
typename TypeTag>
176 return well_ecl_->getAllowCrossFlow();
183 template<
typename TypeTag>
184 const std::vector<bool>&
185 WellInterface<TypeTag>::
197 template<
typename TypeTag>
199 WellInterface<TypeTag>::
200 setWellEfficiencyFactor(
const double efficiency_factor)
202 well_efficiency_factor_ = efficiency_factor;
209 template<
typename TypeTag>
211 WellInterface<TypeTag>::
214 assert(phase_usage_);
216 return *phase_usage_;
223 template<
typename TypeTag>
225 WellInterface<TypeTag>::
226 flowPhaseToEbosCompIdx(
const int phaseIdx )
const
228 const auto& pu = phaseUsage();
229 if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
230 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
231 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
232 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
233 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
234 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
243 template<
typename TypeTag>
245 WellInterface<TypeTag>::
246 flowPhaseToEbosPhaseIdx(
const int phaseIdx )
const
248 const auto& pu = phaseUsage();
249 if (active()[Water] && pu.phase_pos[Water] == phaseIdx) {
250 return FluidSystem::waterPhaseIdx;
252 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) {
253 return FluidSystem::oilPhaseIdx;
255 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) {
256 return FluidSystem::gasPhaseIdx;
259 assert(phaseIdx < 3);
267 template<
typename TypeTag>
269 WellInterface<TypeTag>::
270 numComponents()
const
273 if (number_of_phases_ == 2) {
277 int numComp = FluidSystem::numComponents;
288 template<
typename TypeTag>
290 WellInterface<TypeTag>::
297 WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
298 if (injection.injectorType == WellInjector::GAS) {
299 double solvent_fraction = well_ecl_->getSolventFraction(current_step_);
300 return solvent_fraction;
311 template<
typename TypeTag>
313 WellInterface<TypeTag>::
320 WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
321 WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_);
323 if (injection.injectorType == WellInjector::WATER) {
324 const double polymer_injection_concentration = polymer.m_polymerConcentration;
325 return polymer_injection_concentration;
336 template<
typename TypeTag>
338 WellInterface<TypeTag>::
339 mostStrictBhpFromBhpLimits()
const
344 switch( well_type_ ) {
346 bhp = std::numeric_limits<double>::max();
349 bhp = -std::numeric_limits<double>::max();
352 OPM_THROW(std::logic_error,
"Expected PRODUCER or INJECTOR type for well " << name());
356 const int nwc = well_controls_get_num(well_controls_);
358 for (
int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
360 if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
362 const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
366 if (bhp_target < bhp) {
371 if (bhp_target > bhp) {
376 OPM_THROW(std::logic_error,
"Expected PRODUCER or INJECTOR type for well " << name());
387 template<
typename TypeTag>
389 WellInterface<TypeTag>::
390 wellHasTHPConstraints()
const
392 const int nwc = well_controls_get_num(well_controls_);
393 for (
int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
394 if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
405 template<
typename TypeTag>
407 WellInterface<TypeTag>::
408 updateWellControl(WellState& well_state,
409 wellhelpers::WellSwitchingLogger& logger)
const
411 const int np = number_of_phases_;
412 const int w = index_of_well_;
414 const int old_control_index = well_state.currentControls()[w];
418 WellControls* wc = well_controls_;
423 const int nwc = well_controls_get_num(wc);
425 int current = well_state.currentControls()[w];
427 for (; ctrl_index < nwc; ++ctrl_index) {
428 if (ctrl_index == current) {
434 if (wellhelpers::constraintBroken(
435 well_state.bhp(), well_state.thp(), well_state.wellRates(),
436 w, np, well_type_, wc, ctrl_index)) {
442 if (ctrl_index != nwc) {
444 well_state.currentControls()[w] = ctrl_index;
445 current = well_state.currentControls()[w];
446 well_controls_set_current( wc, current);
450 const int updated_control_index = well_state.currentControls()[w];
453 if (updated_control_index != old_control_index) {
454 logger.wellSwitched(name(),
455 well_controls_iget_type(wc, old_control_index),
456 well_controls_iget_type(wc, updated_control_index));
459 if (updated_control_index != old_control_index) {
460 updateWellStateWithTarget(updated_control_index, well_state);
468 template<
typename TypeTag>
470 WellInterface<TypeTag>::
471 checkRateEconLimits(
const WellEconProductionLimits& econ_production_limits,
472 const WellState& well_state)
const
474 const Opm::PhaseUsage& pu = phaseUsage();
475 const int np = number_of_phases_;
477 if (econ_production_limits.onMinOilRate()) {
478 assert(active()[Oil]);
479 const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
480 const double min_oil_rate = econ_production_limits.minOilRate();
481 if (std::abs(oil_rate) < min_oil_rate) {
486 if (econ_production_limits.onMinGasRate() ) {
487 assert(active()[Gas]);
488 const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
489 const double min_gas_rate = econ_production_limits.minGasRate();
490 if (std::abs(gas_rate) < min_gas_rate) {
495 if (econ_production_limits.onMinLiquidRate() ) {
496 assert(active()[Oil]);
497 assert(active()[Water]);
498 const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
499 const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
500 const double liquid_rate = oil_rate + water_rate;
501 const double min_liquid_rate = econ_production_limits.minLiquidRate();
502 if (std::abs(liquid_rate) < min_liquid_rate) {
507 if (econ_production_limits.onMinReservoirFluidRate()) {
508 OpmLog::warning(
"NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE",
"Minimum reservoir fluid production rate limit is not supported yet");
519 template<
typename TypeTag>
520 typename WellInterface<TypeTag>::RatioCheckTuple
521 WellInterface<TypeTag>::
522 checkMaxWaterCutLimit(
const WellEconProductionLimits& econ_production_limits,
523 const WellState& well_state)
const
525 bool water_cut_limit_violated =
false;
526 int worst_offending_connection = INVALIDCONNECTION;
527 bool last_connection =
false;
528 double violation_extent = -1.0;
530 const int np = number_of_phases_;
531 const Opm::PhaseUsage& pu = phaseUsage();
532 const int well_number = index_of_well_;
534 assert(active()[Oil]);
535 assert(active()[Water]);
537 const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
538 const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
539 const double liquid_rate = oil_rate + water_rate;
541 if (std::abs(liquid_rate) != 0.) {
542 water_cut = water_rate / liquid_rate;
547 const double max_water_cut_limit = econ_production_limits.maxWaterCut();
548 if (water_cut > max_water_cut_limit) {
549 water_cut_limit_violated =
true;
552 if (water_cut_limit_violated) {
554 const int perf_start = first_perf_;
555 const int perf_number = number_of_perforations_;
557 std::vector<double> water_cut_perf(perf_number);
558 for (
int perf = 0; perf < perf_number; ++perf) {
559 const int i_perf = perf_start + perf;
560 const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
561 const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
562 const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
563 if (std::abs(liquid_perf_rate) != 0.) {
564 water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
566 water_cut_perf[perf] = 0.;
570 last_connection = (perf_number == 1);
571 if (last_connection) {
572 worst_offending_connection = 0;
573 violation_extent = water_cut_perf[0] / max_water_cut_limit;
574 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
577 double max_water_cut_perf = 0.;
578 for (
int perf = 0; perf < perf_number; ++perf) {
579 if (water_cut_perf[perf] > max_water_cut_perf) {
580 worst_offending_connection = perf;
581 max_water_cut_perf = water_cut_perf[perf];
585 assert(max_water_cut_perf != 0.);
586 assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
588 violation_extent = max_water_cut_perf / max_water_cut_limit;
591 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
598 template<
typename TypeTag>
599 typename WellInterface<TypeTag>::RatioCheckTuple
600 WellInterface<TypeTag>::
601 checkRatioEconLimits(
const WellEconProductionLimits& econ_production_limits,
602 const WellState& well_state)
const
612 bool any_limit_violated =
false;
613 bool last_connection =
false;
614 int worst_offending_connection = INVALIDCONNECTION;
615 double violation_extent = -1.0;
617 if (econ_production_limits.onMaxWaterCut()) {
618 const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state);
619 bool water_cut_violated = std::get<0>(water_cut_return);
620 if (water_cut_violated) {
621 any_limit_violated =
true;
622 const double violation_extent_water_cut = std::get<3>(water_cut_return);
623 if (violation_extent_water_cut > violation_extent) {
624 violation_extent = violation_extent_water_cut;
625 worst_offending_connection = std::get<2>(water_cut_return);
626 last_connection = std::get<1>(water_cut_return);
631 if (econ_production_limits.onMaxGasOilRatio()) {
632 OpmLog::warning(
"NOT_SUPPORTING_MAX_GOR",
"the support for max Gas-Oil ratio is not implemented yet!");
635 if (econ_production_limits.onMaxWaterGasRatio()) {
636 OpmLog::warning(
"NOT_SUPPORTING_MAX_WGR",
"the support for max Water-Gas ratio is not implemented yet!");
639 if (econ_production_limits.onMaxGasLiquidRatio()) {
640 OpmLog::warning(
"NOT_SUPPORTING_MAX_GLR",
"the support for max Gas-Liquid ratio is not implemented yet!");
643 if (any_limit_violated) {
644 assert(worst_offending_connection >=0);
645 assert(violation_extent > 1.);
648 return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
655 template<
typename TypeTag>
657 WellInterface<TypeTag>::
658 updateListEconLimited(
const WellState& well_state,
659 DynamicListEconLimited& list_econ_limited)
const
662 if (wellType() != PRODUCER) {
667 bool rate_limit_violated =
false;
668 const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_);
671 if ( !econ_production_limits.onAnyEffectiveLimit() ) {
675 const std::string well_name = name();
679 const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
680 if (quantity_limit == WellEcon::POTN) {
681 const std::string msg = std::string(
"POTN limit for well ") + well_name + std::string(
" is not supported for the moment. \n")
682 + std::string(
"All the limits will be evaluated based on RATE. ");
683 OpmLog::warning(
"NOT_SUPPORTING_POTN", msg);
686 if (econ_production_limits.onAnyRateLimit()) {
687 rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state);
690 if (rate_limit_violated) {
691 if (econ_production_limits.endRun()) {
692 const std::string warning_message = std::string(
"ending run after well closed due to economic limits is not supported yet \n")
693 + std::string(
"the program will keep running after ") + well_name + std::string(
" is closed");
694 OpmLog::warning(
"NOT_SUPPORTING_ENDRUN", warning_message);
697 if (econ_production_limits.validFollowonWell()) {
698 OpmLog::warning(
"NOT_SUPPORTING_FOLLOWONWELL",
"opening following on well after well closed is not supported yet");
701 if (well_ecl_->getAutomaticShutIn()) {
702 list_econ_limited.addShutWell(well_name);
703 const std::string msg = std::string(
"well ") + well_name + std::string(
" will be shut in due to rate economic limit");
706 list_econ_limited.addStoppedWell(well_name);
707 const std::string msg = std::string(
"well ") + well_name + std::string(
" will be stopped due to rate economic limit");
715 bool ratio_limits_violated =
false;
716 RatioCheckTuple ratio_check_return;
718 if (econ_production_limits.onAnyRatioLimit()) {
719 ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state);
720 ratio_limits_violated = std::get<0>(ratio_check_return);
723 if (ratio_limits_violated) {
724 const WellEcon::WorkoverEnum workover = econ_production_limits.workover();
728 const bool last_connection = std::get<1>(ratio_check_return);
729 const int worst_offending_connection = std::get<2>(ratio_check_return);
731 assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_));
733 const int cell_worst_offending_connection = well_cells_[worst_offending_connection];
734 list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
735 const std::string msg = std::string(
"Connection ") + std::to_string(worst_offending_connection) + std::string(
" for well ")
736 + well_name + std::string(
" will be closed due to economic limit");
739 if (last_connection) {
741 list_econ_limited.addShutWell(well_name);
742 const std::string msg2 = well_name + std::string(
" will be shut due to the last connection closed");
749 if (well_ecl_->getAutomaticShutIn()) {
750 list_econ_limited.addShutWell(well_name);
751 const std::string msg = well_name + std::string(
" will be shut due to ratio economic limit");
754 list_econ_limited.addStoppedWell(well_name);
755 const std::string msg = well_name + std::string(
" will be stopped due to ratio economic limit");
764 OpmLog::warning(
"NOT_SUPPORTED_WORKOVER_TYPE",
"not supporting workover type " + WellEcon::WorkoverEnumToString(workover) );
774 template<
typename TypeTag>
776 WellInterface<TypeTag>::
777 computeRepRadiusPerfLength(
const Grid& grid,
778 const std::map<int, int>& cartesian_to_compressed)
780 const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
781 auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
782 auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
784 const int nperf = number_of_perforations_;
786 perf_rep_radius_.clear();
787 perf_length_.clear();
788 bore_diameters_.clear();
790 perf_rep_radius_.reserve(nperf);
791 perf_length_.reserve(nperf);
792 bore_diameters_.reserve(nperf);
795 const auto& completionSet = well_ecl_->getCompletions(current_step_);
796 for (
size_t c=0; c<completionSet.size(); c++) {
797 const auto& completion = completionSet.get(c);
798 if (completion.getState() == WellCompletion::OPEN) {
799 const int i = completion.getI();
800 const int j = completion.getJ();
801 const int k = completion.getK();
803 const int* cpgdim = cart_dims;
804 const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
805 const std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
806 if (cgit == cartesian_to_compressed.end()) {
807 OPM_THROW(std::runtime_error,
"Cell with i,j,k indices " << i <<
' ' << j <<
' '
808 << k <<
" not found in grid (well = " << name() <<
')');
810 const int cell = cgit->second;
813 double radius = 0.5*completion.getDiameter();
815 radius = 0.5*unit::feet;
816 OPM_MESSAGE(
"**** Warning: Well bore internal radius set to " << radius);
819 const std::array<double, 3> cubical =
820 WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
822 WellCompletion::DirectionEnum direction = completion.getDirection();
828 case Opm::WellCompletion::DirectionEnum::X:
829 re = std::sqrt(cubical[1] * cubical[2] / M_PI);
830 perf_length = cubical[0];
832 case Opm::WellCompletion::DirectionEnum::Y:
833 re = std::sqrt(cubical[0] * cubical[2] / M_PI);
834 perf_length = cubical[1];
836 case Opm::WellCompletion::DirectionEnum::Z:
837 re = std::sqrt(cubical[0] * cubical[1] / M_PI);
838 perf_length = cubical[2];
841 OPM_THROW(std::runtime_error,
" Dirtecion of well is not supported ");
844 const double repR = std::sqrt(re * radius);
845 perf_rep_radius_.push_back(repR);
846 perf_length_.push_back(perf_length);
847 bore_diameters_.push_back(2. * radius);
const std::string & name() const
Well name.
Definition: WellInterface_impl.hpp:138
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
Definition: WellInterface.hpp:60
WellControls * wellControls() const
Well controls.
Definition: WellInterface_impl.hpp:162
WellInterface(const Well *well, const int time_step, const Wells *wells, const ModelParameters ¶m)
Constructor.
Definition: WellInterface_impl.hpp:28
WellType wellType() const
Well type, INJECTOR or PRODUCER.
Definition: WellInterface_impl.hpp:150