00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 namespace Opm
00023 {
00024
00025
00026 template<typename TypeTag>
00027 WellInterface<TypeTag>::
00028 WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
00029 : well_ecl_(well)
00030 , current_step_(time_step)
00031 , param_(param)
00032 {
00033 if (!well) {
00034 OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface");
00035 }
00036
00037 if (time_step < 0) {
00038 OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
00039 }
00040
00041 if (!wells) {
00042 OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface");
00043 }
00044
00045 const std::string& well_name = well->name();
00046
00047
00048 int index_well;
00049 for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
00050 if (well_name == std::string(wells->name[index_well])) {
00051 break;
00052 }
00053 }
00054
00055
00056
00057 assert(index_well != wells->number_of_wells);
00058
00059 index_of_well_ = index_well;
00060 well_type_ = wells->type[index_well];
00061 number_of_phases_ = wells->number_of_phases;
00062
00063
00064 {
00065 comp_frac_.resize(number_of_phases_);
00066 const int index_begin = index_well * number_of_phases_;
00067 std::copy(wells->comp_frac + index_begin,
00068 wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
00069 }
00070
00071 well_controls_ = wells->ctrls[index_well];
00072
00073 ref_depth_ = wells->depth_ref[index_well];
00074
00075
00076 {
00077 const int perf_index_begin = wells->well_connpos[index_well];
00078 const int perf_index_end = wells->well_connpos[index_well + 1];
00079 number_of_perforations_ = perf_index_end - perf_index_begin;
00080 first_perf_ = perf_index_begin;
00081
00082 well_cells_.resize(number_of_perforations_);
00083 std::copy(wells->well_cells + perf_index_begin,
00084 wells->well_cells + perf_index_end,
00085 well_cells_.begin() );
00086
00087 well_index_.resize(number_of_perforations_);
00088 std::copy(wells->WI + perf_index_begin,
00089 wells->WI + perf_index_end,
00090 well_index_.begin() );
00091
00092 saturation_table_number_.resize(number_of_perforations_);
00093 std::copy(wells->sat_table_id + perf_index_begin,
00094 wells->sat_table_id + perf_index_end,
00095 saturation_table_number_.begin() );
00096 }
00097
00098 well_efficiency_factor_ = 1.0;
00099 }
00100
00101
00102
00103
00104
00105 template<typename TypeTag>
00106 void
00107 WellInterface<TypeTag>::
00108 init(const PhaseUsage* phase_usage_arg,
00109 const std::vector<bool>* active_arg,
00110 const std::vector<double>& ,
00111 const double gravity_arg,
00112 const int )
00113 {
00114 phase_usage_ = phase_usage_arg;
00115 active_ = active_arg;
00116 gravity_ = gravity_arg;
00117 }
00118
00119
00120
00121
00122
00123 template<typename TypeTag>
00124 void
00125 WellInterface<TypeTag>::
00126 setVFPProperties(const VFPProperties* vfp_properties_arg)
00127 {
00128 vfp_properties_ = vfp_properties_arg;
00129 }
00130
00131
00132
00133
00134
00135 template<typename TypeTag>
00136 const std::string&
00137 WellInterface<TypeTag>::
00138 name() const
00139 {
00140 return well_ecl_->name();
00141 }
00142
00143
00144
00145
00146
00147 template<typename TypeTag>
00148 WellType
00149 WellInterface<TypeTag>::
00150 wellType() const
00151 {
00152 return well_type_;
00153 }
00154
00155
00156
00157
00158
00159 template<typename TypeTag>
00160 WellControls*
00161 WellInterface<TypeTag>::
00162 wellControls() const
00163 {
00164 return well_controls_;
00165 }
00166
00167
00168
00169
00170
00171 template<typename TypeTag>
00172 bool
00173 WellInterface<TypeTag>::
00174 getAllowCrossFlow() const
00175 {
00176 return well_ecl_->getAllowCrossFlow();
00177 }
00178
00179
00180
00181
00182
00183 template<typename TypeTag>
00184 const std::vector<bool>&
00185 WellInterface<TypeTag>::
00186 active() const
00187 {
00188 assert(active_);
00189
00190 return *active_;
00191 }
00192
00193
00194
00195
00196
00197 template<typename TypeTag>
00198 void
00199 WellInterface<TypeTag>::
00200 setWellEfficiencyFactor(const double efficiency_factor)
00201 {
00202 well_efficiency_factor_ = efficiency_factor;
00203 }
00204
00205
00206
00207
00208
00209 template<typename TypeTag>
00210 const PhaseUsage&
00211 WellInterface<TypeTag>::
00212 phaseUsage() const
00213 {
00214 assert(phase_usage_);
00215
00216 return *phase_usage_;
00217 }
00218
00219
00220
00221
00222
00223 template<typename TypeTag>
00224 int
00225 WellInterface<TypeTag>::
00226 flowPhaseToEbosCompIdx( const int phaseIdx ) const
00227 {
00228 const auto& pu = phaseUsage();
00229 if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
00230 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
00231 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
00232 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
00233 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
00234 return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
00235
00236
00237 return phaseIdx;
00238 }
00239
00240
00241
00242
00243 template<typename TypeTag>
00244 int
00245 WellInterface<TypeTag>::
00246 flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
00247 {
00248 const auto& pu = phaseUsage();
00249 if (active()[Water] && pu.phase_pos[Water] == phaseIdx) {
00250 return FluidSystem::waterPhaseIdx;
00251 }
00252 if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) {
00253 return FluidSystem::oilPhaseIdx;
00254 }
00255 if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) {
00256 return FluidSystem::gasPhaseIdx;
00257 }
00258
00259 assert(phaseIdx < 3);
00260
00261 return phaseIdx;
00262 }
00263
00264
00265
00266
00267 template<typename TypeTag>
00268 int
00269 WellInterface<TypeTag>::
00270 numComponents() const
00271 {
00272
00273 if (number_of_phases_ == 2) {
00274 return 2;
00275 }
00276
00277 int numComp = FluidSystem::numComponents;
00278
00279 if (has_solvent) {
00280 numComp++;
00281 }
00282 return numComp;
00283 }
00284
00285
00286
00287
00288 template<typename TypeTag>
00289 double
00290 WellInterface<TypeTag>::
00291 wsolvent() const
00292 {
00293 if (!has_solvent) {
00294 return 0.0;
00295 }
00296
00297 WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
00298 if (injection.injectorType == WellInjector::GAS) {
00299 double solvent_fraction = well_ecl_->getSolventFraction(current_step_);
00300 return solvent_fraction;
00301 }
00302
00303 assert(false);
00304 return 0.0;
00305 }
00306
00307
00308
00309
00310
00311 template<typename TypeTag>
00312 double
00313 WellInterface<TypeTag>::
00314 wpolymer() const
00315 {
00316 if (!has_polymer) {
00317 return 0.0;
00318 }
00319
00320 WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
00321 WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_);
00322
00323 if (injection.injectorType == WellInjector::WATER) {
00324 const double polymer_injection_concentration = polymer.m_polymerConcentration;
00325 return polymer_injection_concentration;
00326 }
00327
00328 assert(false);
00329 return 0.0;
00330 }
00331
00332
00333
00334
00335
00336 template<typename TypeTag>
00337 double
00338 WellInterface<TypeTag>::
00339 mostStrictBhpFromBhpLimits() const
00340 {
00341 double bhp;
00342
00343
00344 switch( well_type_ ) {
00345 case INJECTOR:
00346 bhp = std::numeric_limits<double>::max();
00347 break;
00348 case PRODUCER:
00349 bhp = -std::numeric_limits<double>::max();
00350 break;
00351 default:
00352 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
00353 }
00354
00355
00356 const int nwc = well_controls_get_num(well_controls_);
00357
00358 for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
00359
00360 if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
00361
00362 const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
00363
00364 switch(well_type_) {
00365 case INJECTOR:
00366 if (bhp_target < bhp) {
00367 bhp = bhp_target;
00368 }
00369 break;
00370 case PRODUCER:
00371 if (bhp_target > bhp) {
00372 bhp = bhp_target;
00373 }
00374 break;
00375 default:
00376 OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
00377 }
00378 }
00379 }
00380
00381 return bhp;
00382 }
00383
00384
00385
00386
00387 template<typename TypeTag>
00388 bool
00389 WellInterface<TypeTag>::
00390 wellHasTHPConstraints() const
00391 {
00392 const int nwc = well_controls_get_num(well_controls_);
00393 for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
00394 if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
00395 return true;
00396 }
00397 }
00398 return false;
00399 }
00400
00401
00402
00403
00404
00405 template<typename TypeTag>
00406 void
00407 WellInterface<TypeTag>::
00408 updateWellControl(WellState& well_state,
00409 wellhelpers::WellSwitchingLogger& logger) const
00410 {
00411 const int np = number_of_phases_;
00412 const int w = index_of_well_;
00413
00414 const int old_control_index = well_state.currentControls()[w];
00415
00416
00417
00418 WellControls* wc = well_controls_;
00419
00420
00421
00422
00423 const int nwc = well_controls_get_num(wc);
00424
00425 int current = well_state.currentControls()[w];
00426 int ctrl_index = 0;
00427 for (; ctrl_index < nwc; ++ctrl_index) {
00428 if (ctrl_index == current) {
00429
00430
00431
00432 continue;
00433 }
00434 if (wellhelpers::constraintBroken(
00435 well_state.bhp(), well_state.thp(), well_state.wellRates(),
00436 w, np, well_type_, wc, ctrl_index)) {
00437
00438 break;
00439 }
00440 }
00441
00442 if (ctrl_index != nwc) {
00443
00444 well_state.currentControls()[w] = ctrl_index;
00445 current = well_state.currentControls()[w];
00446 well_controls_set_current( wc, current);
00447 }
00448
00449
00450 const int updated_control_index = well_state.currentControls()[w];
00451
00452
00453 if (updated_control_index != old_control_index) {
00454 logger.wellSwitched(name(),
00455 well_controls_iget_type(wc, old_control_index),
00456 well_controls_iget_type(wc, updated_control_index));
00457 }
00458
00459 if (updated_control_index != old_control_index) {
00460 updateWellStateWithTarget(updated_control_index, well_state);
00461 }
00462 }
00463
00464
00465
00466
00467
00468 template<typename TypeTag>
00469 bool
00470 WellInterface<TypeTag>::
00471 checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
00472 const WellState& well_state) const
00473 {
00474 const Opm::PhaseUsage& pu = phaseUsage();
00475 const int np = number_of_phases_;
00476
00477 if (econ_production_limits.onMinOilRate()) {
00478 assert(active()[Oil]);
00479 const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
00480 const double min_oil_rate = econ_production_limits.minOilRate();
00481 if (std::abs(oil_rate) < min_oil_rate) {
00482 return true;
00483 }
00484 }
00485
00486 if (econ_production_limits.onMinGasRate() ) {
00487 assert(active()[Gas]);
00488 const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
00489 const double min_gas_rate = econ_production_limits.minGasRate();
00490 if (std::abs(gas_rate) < min_gas_rate) {
00491 return true;
00492 }
00493 }
00494
00495 if (econ_production_limits.onMinLiquidRate() ) {
00496 assert(active()[Oil]);
00497 assert(active()[Water]);
00498 const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
00499 const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
00500 const double liquid_rate = oil_rate + water_rate;
00501 const double min_liquid_rate = econ_production_limits.minLiquidRate();
00502 if (std::abs(liquid_rate) < min_liquid_rate) {
00503 return true;
00504 }
00505 }
00506
00507 if (econ_production_limits.onMinReservoirFluidRate()) {
00508 OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
00509 }
00510
00511 return false;
00512 }
00513
00514
00515
00516
00517
00518
00519 template<typename TypeTag>
00520 typename WellInterface<TypeTag>::RatioCheckTuple
00521 WellInterface<TypeTag>::
00522 checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
00523 const WellState& well_state) const
00524 {
00525 bool water_cut_limit_violated = false;
00526 int worst_offending_connection = INVALIDCONNECTION;
00527 bool last_connection = false;
00528 double violation_extent = -1.0;
00529
00530 const int np = number_of_phases_;
00531 const Opm::PhaseUsage& pu = phaseUsage();
00532 const int well_number = index_of_well_;
00533
00534 assert(active()[Oil]);
00535 assert(active()[Water]);
00536
00537 const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
00538 const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
00539 const double liquid_rate = oil_rate + water_rate;
00540 double water_cut;
00541 if (std::abs(liquid_rate) != 0.) {
00542 water_cut = water_rate / liquid_rate;
00543 } else {
00544 water_cut = 0.0;
00545 }
00546
00547 const double max_water_cut_limit = econ_production_limits.maxWaterCut();
00548 if (water_cut > max_water_cut_limit) {
00549 water_cut_limit_violated = true;
00550 }
00551
00552 if (water_cut_limit_violated) {
00553
00554 const int perf_start = first_perf_;
00555 const int perf_number = number_of_perforations_;
00556
00557 std::vector<double> water_cut_perf(perf_number);
00558 for (int perf = 0; perf < perf_number; ++perf) {
00559 const int i_perf = perf_start + perf;
00560 const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
00561 const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
00562 const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
00563 if (std::abs(liquid_perf_rate) != 0.) {
00564 water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
00565 } else {
00566 water_cut_perf[perf] = 0.;
00567 }
00568 }
00569
00570 last_connection = (perf_number == 1);
00571 if (last_connection) {
00572 worst_offending_connection = 0;
00573 violation_extent = water_cut_perf[0] / max_water_cut_limit;
00574 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
00575 }
00576
00577 double max_water_cut_perf = 0.;
00578 for (int perf = 0; perf < perf_number; ++perf) {
00579 if (water_cut_perf[perf] > max_water_cut_perf) {
00580 worst_offending_connection = perf;
00581 max_water_cut_perf = water_cut_perf[perf];
00582 }
00583 }
00584
00585 assert(max_water_cut_perf != 0.);
00586 assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
00587
00588 violation_extent = max_water_cut_perf / max_water_cut_limit;
00589 }
00590
00591 return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
00592 }
00593
00594
00595
00596
00597
00598 template<typename TypeTag>
00599 typename WellInterface<TypeTag>::RatioCheckTuple
00600 WellInterface<TypeTag>::
00601 checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
00602 const WellState& well_state) const
00603 {
00604
00605
00606
00607
00608
00609
00610
00611
00612 bool any_limit_violated = false;
00613 bool last_connection = false;
00614 int worst_offending_connection = INVALIDCONNECTION;
00615 double violation_extent = -1.0;
00616
00617 if (econ_production_limits.onMaxWaterCut()) {
00618 const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state);
00619 bool water_cut_violated = std::get<0>(water_cut_return);
00620 if (water_cut_violated) {
00621 any_limit_violated = true;
00622 const double violation_extent_water_cut = std::get<3>(water_cut_return);
00623 if (violation_extent_water_cut > violation_extent) {
00624 violation_extent = violation_extent_water_cut;
00625 worst_offending_connection = std::get<2>(water_cut_return);
00626 last_connection = std::get<1>(water_cut_return);
00627 }
00628 }
00629 }
00630
00631 if (econ_production_limits.onMaxGasOilRatio()) {
00632 OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
00633 }
00634
00635 if (econ_production_limits.onMaxWaterGasRatio()) {
00636 OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
00637 }
00638
00639 if (econ_production_limits.onMaxGasLiquidRatio()) {
00640 OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
00641 }
00642
00643 if (any_limit_violated) {
00644 assert(worst_offending_connection >=0);
00645 assert(violation_extent > 1.);
00646 }
00647
00648 return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
00649 }
00650
00651
00652
00653
00654
00655 template<typename TypeTag>
00656 void
00657 WellInterface<TypeTag>::
00658 updateListEconLimited(const WellState& well_state,
00659 DynamicListEconLimited& list_econ_limited) const
00660 {
00661
00662 if (wellType() != PRODUCER) {
00663 return;
00664 }
00665
00666
00667 bool rate_limit_violated = false;
00668 const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_);
00669
00670
00671 if ( !econ_production_limits.onAnyEffectiveLimit() ) {
00672 return;
00673 }
00674
00675 const std::string well_name = name();
00676
00677
00678
00679 const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
00680 if (quantity_limit == WellEcon::POTN) {
00681 const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n")
00682 + std::string("All the limits will be evaluated based on RATE. ");
00683 OpmLog::warning("NOT_SUPPORTING_POTN", msg);
00684 }
00685
00686 if (econ_production_limits.onAnyRateLimit()) {
00687 rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state);
00688 }
00689
00690 if (rate_limit_violated) {
00691 if (econ_production_limits.endRun()) {
00692 const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n")
00693 + std::string("the program will keep running after ") + well_name + std::string(" is closed");
00694 OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message);
00695 }
00696
00697 if (econ_production_limits.validFollowonWell()) {
00698 OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
00699 }
00700
00701 if (well_ecl_->getAutomaticShutIn()) {
00702 list_econ_limited.addShutWell(well_name);
00703 const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to rate economic limit");
00704 OpmLog::info(msg);
00705 } else {
00706 list_econ_limited.addStoppedWell(well_name);
00707 const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to rate economic limit");
00708 OpmLog::info(msg);
00709 }
00710
00711 return;
00712 }
00713
00714
00715 bool ratio_limits_violated = false;
00716 RatioCheckTuple ratio_check_return;
00717
00718 if (econ_production_limits.onAnyRatioLimit()) {
00719 ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state);
00720 ratio_limits_violated = std::get<0>(ratio_check_return);
00721 }
00722
00723 if (ratio_limits_violated) {
00724 const WellEcon::WorkoverEnum workover = econ_production_limits.workover();
00725 switch (workover) {
00726 case WellEcon::CON:
00727 {
00728 const bool last_connection = std::get<1>(ratio_check_return);
00729 const int worst_offending_connection = std::get<2>(ratio_check_return);
00730
00731 assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_));
00732
00733 const int cell_worst_offending_connection = well_cells_[worst_offending_connection];
00734 list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
00735 const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ")
00736 + well_name + std::string(" will be closed due to economic limit");
00737 OpmLog::info(msg);
00738
00739 if (last_connection) {
00740
00741 list_econ_limited.addShutWell(well_name);
00742 const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed");
00743 OpmLog::info(msg2);
00744 }
00745 break;
00746 }
00747 case WellEcon::WELL:
00748 {
00749 if (well_ecl_->getAutomaticShutIn()) {
00750 list_econ_limited.addShutWell(well_name);
00751 const std::string msg = well_name + std::string(" will be shut due to ratio economic limit");
00752 OpmLog::info(msg);
00753 } else {
00754 list_econ_limited.addStoppedWell(well_name);
00755 const std::string msg = well_name + std::string(" will be stopped due to ratio economic limit");
00756 OpmLog::info(msg);
00757 }
00758 break;
00759 }
00760 case WellEcon::NONE:
00761 break;
00762 default:
00763 {
00764 OpmLog::warning("NOT_SUPPORTED_WORKOVER_TYPE", "not supporting workover type " + WellEcon::WorkoverEnumToString(workover) );
00765 }
00766 }
00767 }
00768 }
00769
00770
00771
00772
00773
00774 template<typename TypeTag>
00775 void
00776 WellInterface<TypeTag>::
00777 computeRepRadiusPerfLength(const Grid& grid,
00778 const std::map<int, int>& cartesian_to_compressed)
00779 {
00780 const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
00781 auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
00782 auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
00783
00784 const int nperf = number_of_perforations_;
00785
00786 perf_rep_radius_.clear();
00787 perf_length_.clear();
00788 bore_diameters_.clear();
00789
00790 perf_rep_radius_.reserve(nperf);
00791 perf_length_.reserve(nperf);
00792 bore_diameters_.reserve(nperf);
00793
00794
00795 const auto& completionSet = well_ecl_->getCompletions(current_step_);
00796 for (size_t c=0; c<completionSet.size(); c++) {
00797 const auto& completion = completionSet.get(c);
00798 if (completion.getState() == WellCompletion::OPEN) {
00799 const int i = completion.getI();
00800 const int j = completion.getJ();
00801 const int k = completion.getK();
00802
00803 const int* cpgdim = cart_dims;
00804 const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
00805 const std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
00806 if (cgit == cartesian_to_compressed.end()) {
00807 OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
00808 << k << " not found in grid (well = " << name() << ')');
00809 }
00810 const int cell = cgit->second;
00811
00812 {
00813 double radius = 0.5*completion.getDiameter();
00814 if (radius <= 0.0) {
00815 radius = 0.5*unit::feet;
00816 OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
00817 }
00818
00819 const std::array<double, 3> cubical =
00820 WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
00821
00822 WellCompletion::DirectionEnum direction = completion.getDirection();
00823
00824 double re;
00825 double perf_length;
00826
00827 switch (direction) {
00828 case Opm::WellCompletion::DirectionEnum::X:
00829 re = std::sqrt(cubical[1] * cubical[2] / M_PI);
00830 perf_length = cubical[0];
00831 break;
00832 case Opm::WellCompletion::DirectionEnum::Y:
00833 re = std::sqrt(cubical[0] * cubical[2] / M_PI);
00834 perf_length = cubical[1];
00835 break;
00836 case Opm::WellCompletion::DirectionEnum::Z:
00837 re = std::sqrt(cubical[0] * cubical[1] / M_PI);
00838 perf_length = cubical[2];
00839 break;
00840 default:
00841 OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
00842 }
00843
00844 const double repR = std::sqrt(re * radius);
00845 perf_rep_radius_.push_back(repR);
00846 perf_length_.push_back(perf_length);
00847 bore_diameters_.push_back(2. * radius);
00848 }
00849 }
00850 }
00851 }
00852
00853 }