All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
WellInterface_impl.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 namespace Opm
23 {
24 
25 
26  template<typename TypeTag>
28  WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
29  : well_ecl_(well)
30  , current_step_(time_step)
31  , param_(param)
32  {
33  if (!well) {
34  OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface");
35  }
36 
37  if (time_step < 0) {
38  OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
39  }
40 
41  if (!wells) {
42  OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface");
43  }
44 
45  const std::string& well_name = well->name();
46 
47  // looking for the location of the well in the wells struct
48  int index_well;
49  for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
50  if (well_name == std::string(wells->name[index_well])) {
51  break;
52  }
53  }
54 
55  // should not enter the constructor if the well does not exist in the wells struct
56  // here, just another assertion.
57  assert(index_well != wells->number_of_wells);
58 
59  index_of_well_ = index_well;
60  well_type_ = wells->type[index_well];
61  number_of_phases_ = wells->number_of_phases;
62 
63  // copying the comp_frac
64  {
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() );
69  }
70 
71  well_controls_ = wells->ctrls[index_well];
72 
73  ref_depth_ = wells->depth_ref[index_well];
74 
75  // perforations related
76  {
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;
81 
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() );
86 
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() );
91 
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() );
96  }
97 
98  well_efficiency_factor_ = 1.0;
99  }
100 
101 
102 
103 
104 
105  template<typename TypeTag>
106  void
108  init(const PhaseUsage* phase_usage_arg,
109  const std::vector<bool>* active_arg,
110  const std::vector<double>& /* depth_arg */,
111  const double gravity_arg,
112  const int /* num_cells */)
113  {
114  phase_usage_ = phase_usage_arg;
115  active_ = active_arg;
116  gravity_ = gravity_arg;
117  }
118 
119 
120 
121 
122 
123  template<typename TypeTag>
124  void
125  WellInterface<TypeTag>::
126  setVFPProperties(const VFPProperties* vfp_properties_arg)
127  {
128  vfp_properties_ = vfp_properties_arg;
129  }
130 
131 
132 
133 
134 
135  template<typename TypeTag>
136  const std::string&
138  name() const
139  {
140  return well_ecl_->name();
141  }
142 
143 
144 
145 
146 
147  template<typename TypeTag>
148  WellType
150  wellType() const
151  {
152  return well_type_;
153  }
154 
155 
156 
157 
158 
159  template<typename TypeTag>
160  WellControls*
163  {
164  return well_controls_;
165  }
166 
167 
168 
169 
170 
171  template<typename TypeTag>
172  bool
174  getAllowCrossFlow() const
175  {
176  return well_ecl_->getAllowCrossFlow();
177  }
178 
179 
180 
181 
182 
183  template<typename TypeTag>
184  const std::vector<bool>&
185  WellInterface<TypeTag>::
186  active() const
187  {
188  assert(active_);
189 
190  return *active_;
191  }
192 
193 
194 
195 
196 
197  template<typename TypeTag>
198  void
199  WellInterface<TypeTag>::
200  setWellEfficiencyFactor(const double efficiency_factor)
201  {
202  well_efficiency_factor_ = efficiency_factor;
203  }
204 
205 
206 
207 
208 
209  template<typename TypeTag>
210  const PhaseUsage&
211  WellInterface<TypeTag>::
212  phaseUsage() const
213  {
214  assert(phase_usage_);
215 
216  return *phase_usage_;
217  }
218 
219 
220 
221 
222 
223  template<typename TypeTag>
224  int
225  WellInterface<TypeTag>::
226  flowPhaseToEbosCompIdx( const int phaseIdx ) const
227  {
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);
235 
236  // for other phases return the index
237  return phaseIdx;
238  }
239 
240 
241 
242 
243  template<typename TypeTag>
244  int
245  WellInterface<TypeTag>::
246  flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
247  {
248  const auto& pu = phaseUsage();
249  if (active()[Water] && pu.phase_pos[Water] == phaseIdx) {
250  return FluidSystem::waterPhaseIdx;
251  }
252  if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) {
253  return FluidSystem::oilPhaseIdx;
254  }
255  if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) {
256  return FluidSystem::gasPhaseIdx;
257  }
258 
259  assert(phaseIdx < 3);
260  // for other phases return the index
261  return phaseIdx;
262  }
263 
264 
265 
266 
267  template<typename TypeTag>
268  int
269  WellInterface<TypeTag>::
270  numComponents() const
271  {
272  // TODO: how about two phase polymer
273  if (number_of_phases_ == 2) {
274  return 2;
275  }
276 
277  int numComp = FluidSystem::numComponents;
278 
279  if (has_solvent) {
280  numComp++;
281  }
282  return numComp;
283  }
284 
285 
286 
287 
288  template<typename TypeTag>
289  double
290  WellInterface<TypeTag>::
291  wsolvent() const
292  {
293  if (!has_solvent) {
294  return 0.0;
295  }
296 
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;
301  }
302 
303  assert(false);
304  return 0.0;
305  }
306 
307 
308 
309 
310 
311  template<typename TypeTag>
312  double
313  WellInterface<TypeTag>::
314  wpolymer() const
315  {
316  if (!has_polymer) {
317  return 0.0;
318  }
319 
320  WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
321  WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_);
322 
323  if (injection.injectorType == WellInjector::WATER) {
324  const double polymer_injection_concentration = polymer.m_polymerConcentration;
325  return polymer_injection_concentration;
326  }
327 
328  assert(false); // TODO: find a more logical way to handle this situation
329  return 0.0;
330  }
331 
332 
333 
334 
335 
336  template<typename TypeTag>
337  double
338  WellInterface<TypeTag>::
339  mostStrictBhpFromBhpLimits() const
340  {
341  double bhp;
342 
343  // initial bhp value, making the value not usable
344  switch( well_type_ ) {
345  case INJECTOR:
346  bhp = std::numeric_limits<double>::max();
347  break;
348  case PRODUCER:
349  bhp = -std::numeric_limits<double>::max();
350  break;
351  default:
352  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
353  }
354 
355  // The number of the well controls/constraints
356  const int nwc = well_controls_get_num(well_controls_);
357 
358  for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
359  // finding a BHP constraint
360  if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
361  // get the bhp constraint value, it should always be postive assummingly
362  const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
363 
364  switch(well_type_) {
365  case INJECTOR: // using the lower bhp contraint from Injectors
366  if (bhp_target < bhp) {
367  bhp = bhp_target;
368  }
369  break;
370  case PRODUCER:
371  if (bhp_target > bhp) {
372  bhp = bhp_target;
373  }
374  break;
375  default:
376  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
377  } // end of switch
378  }
379  }
380 
381  return bhp;
382  }
383 
384 
385 
386 
387  template<typename TypeTag>
388  bool
389  WellInterface<TypeTag>::
390  wellHasTHPConstraints() const
391  {
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) {
395  return true;
396  }
397  }
398  return false;
399  }
400 
401 
402 
403 
404 
405  template<typename TypeTag>
406  void
407  WellInterface<TypeTag>::
408  updateWellControl(WellState& well_state,
409  wellhelpers::WellSwitchingLogger& logger) const
410  {
411  const int np = number_of_phases_;
412  const int w = index_of_well_;
413 
414  const int old_control_index = well_state.currentControls()[w];
415 
416  // Find, for each well, if any constraints are broken. If so,
417  // switch control to first broken constraint.
418  WellControls* wc = well_controls_;
419 
420  // Loop over all controls except the current one, and also
421  // skip any RESERVOIR_RATE controls, since we cannot
422  // handle those.
423  const int nwc = well_controls_get_num(wc);
424  // the current control index
425  int current = well_state.currentControls()[w];
426  int ctrl_index = 0;
427  for (; ctrl_index < nwc; ++ctrl_index) {
428  if (ctrl_index == current) {
429  // This is the currently used control, so it is
430  // used as an equation. So this is not used as an
431  // inequality constraint, and therefore skipped.
432  continue;
433  }
434  if (wellhelpers::constraintBroken(
435  well_state.bhp(), well_state.thp(), well_state.wellRates(),
436  w, np, well_type_, wc, ctrl_index)) {
437  // ctrl_index will be the index of the broken constraint after the loop.
438  break;
439  }
440  }
441 
442  if (ctrl_index != nwc) {
443  // Constraint number ctrl_index was broken, switch to it.
444  well_state.currentControls()[w] = ctrl_index;
445  current = well_state.currentControls()[w];
446  well_controls_set_current( wc, current);
447  }
448 
449  // the new well control indices after all the related updates,
450  const int updated_control_index = well_state.currentControls()[w];
451 
452  // checking whether control changed
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));
457  }
458 
459  if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) {
460  updateWellStateWithTarget(updated_control_index, well_state);
461  }
462  }
463 
464 
465 
466 
467 
468  template<typename TypeTag>
469  bool
470  WellInterface<TypeTag>::
471  checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
472  const WellState& well_state) const
473  {
474  const Opm::PhaseUsage& pu = phaseUsage();
475  const int np = number_of_phases_;
476 
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) {
482  return true;
483  }
484  }
485 
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) {
491  return true;
492  }
493  }
494 
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) {
503  return true;
504  }
505  }
506 
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");
509  }
510 
511  return false;
512  }
513 
514 
515 
516 
517 
518 
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
524  {
525  bool water_cut_limit_violated = false;
526  int worst_offending_connection = INVALIDCONNECTION;
527  bool last_connection = false;
528  double violation_extent = -1.0;
529 
530  const int np = number_of_phases_;
531  const Opm::PhaseUsage& pu = phaseUsage();
532  const int well_number = index_of_well_;
533 
534  assert(active()[Oil]);
535  assert(active()[Water]);
536 
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;
540  double water_cut;
541  if (std::abs(liquid_rate) != 0.) {
542  water_cut = water_rate / liquid_rate;
543  } else {
544  water_cut = 0.0;
545  }
546 
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;
550  }
551 
552  if (water_cut_limit_violated) {
553  // need to handle the worst_offending_connection
554  const int perf_start = first_perf_;
555  const int perf_number = number_of_perforations_;
556 
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;
565  } else {
566  water_cut_perf[perf] = 0.;
567  }
568  }
569 
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);
575  }
576 
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];
582  }
583  }
584 
585  assert(max_water_cut_perf != 0.);
586  assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
587 
588  violation_extent = max_water_cut_perf / max_water_cut_limit;
589  }
590 
591  return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
592  }
593 
594 
595 
596 
597 
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
603  {
604  // TODO: not sure how to define the worst-offending connection when more than one
605  // ratio related limit is violated.
606  // The defintion used here is that we define the violation extent based on the
607  // ratio between the value and the corresponding limit.
608  // For each violated limit, we decide the worst-offending connection separately.
609  // Among the worst-offending connections, we use the one has the biggest violation
610  // extent.
611 
612  bool any_limit_violated = false;
613  bool last_connection = false;
614  int worst_offending_connection = INVALIDCONNECTION;
615  double violation_extent = -1.0;
616 
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);
627  }
628  }
629  }
630 
631  if (econ_production_limits.onMaxGasOilRatio()) {
632  OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
633  }
634 
635  if (econ_production_limits.onMaxWaterGasRatio()) {
636  OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
637  }
638 
639  if (econ_production_limits.onMaxGasLiquidRatio()) {
640  OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
641  }
642 
643  if (any_limit_violated) {
644  assert(worst_offending_connection >=0);
645  assert(violation_extent > 1.);
646  }
647 
648  return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
649  }
650 
651 
652 
653 
654 
655  template<typename TypeTag>
656  void
657  WellInterface<TypeTag>::
658  updateListEconLimited(const WellState& well_state,
659  DynamicListEconLimited& list_econ_limited) const
660  {
661  // economic limits only apply for production wells.
662  if (wellType() != PRODUCER) {
663  return;
664  }
665 
666  // flag to check if the mim oil/gas rate limit is violated
667  bool rate_limit_violated = false;
668  const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_);
669 
670  // if no limit is effective here, then continue to the next well
671  if ( !econ_production_limits.onAnyEffectiveLimit() ) {
672  return;
673  }
674 
675  const std::string well_name = name();
676 
677  // for the moment, we only handle rate limits, not handling potential limits
678  // the potential limits should not be difficult to add
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);
684  }
685 
686  if (econ_production_limits.onAnyRateLimit()) {
687  rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state);
688  }
689 
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);
695  }
696 
697  if (econ_production_limits.validFollowonWell()) {
698  OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
699  }
700 
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");
704  OpmLog::info(msg);
705  } else {
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");
708  OpmLog::info(msg);
709  }
710  // the well is closed, not need to check other limits
711  return;
712  }
713 
714  // checking for ratio related limits, mostly all kinds of ratio.
715  bool ratio_limits_violated = false;
716  RatioCheckTuple ratio_check_return;
717 
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);
721  }
722 
723  if (ratio_limits_violated) {
724  const WellEcon::WorkoverEnum workover = econ_production_limits.workover();
725  switch (workover) {
726  case WellEcon::CON:
727  {
728  const bool last_connection = std::get<1>(ratio_check_return);
729  const int worst_offending_connection = std::get<2>(ratio_check_return);
730 
731  assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_));
732 
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");
737  OpmLog::info(msg);
738 
739  if (last_connection) {
740  // TODO: there is more things to check here
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");
743  OpmLog::info(msg2);
744  }
745  break;
746  }
747  case WellEcon::WELL:
748  {
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");
752  OpmLog::info(msg);
753  } else {
754  list_econ_limited.addStoppedWell(well_name);
755  const std::string msg = well_name + std::string(" will be stopped due to ratio economic limit");
756  OpmLog::info(msg);
757  }
758  break;
759  }
760  case WellEcon::NONE:
761  break;
762  default:
763  {
764  OpmLog::warning("NOT_SUPPORTED_WORKOVER_TYPE", "not supporting workover type " + WellEcon::WorkoverEnumToString(workover) );
765  }
766  }
767  }
768  }
769 
770 
771 
772 
773 
774  template<typename TypeTag>
775  void
776  WellInterface<TypeTag>::
777  computeRepRadiusPerfLength(const Grid& grid,
778  const std::map<int, int>& cartesian_to_compressed)
779  {
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);
783 
784  const int nperf = number_of_perforations_;
785 
786  perf_rep_radius_.clear();
787  perf_length_.clear();
788  bore_diameters_.clear();
789 
790  perf_rep_radius_.reserve(nperf);
791  perf_length_.reserve(nperf);
792  bore_diameters_.reserve(nperf);
793 
794  // COMPDAT handling
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();
802 
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() << ')');
809  }
810  const int cell = cgit->second;
811 
812  {
813  double radius = 0.5*completion.getDiameter();
814  if (radius <= 0.0) {
815  radius = 0.5*unit::feet;
816  OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
817  }
818 
819  const std::array<double, 3> cubical =
820  WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
821 
822  WellCompletion::DirectionEnum direction = completion.getDirection();
823 
824  double re; // area equivalent radius of the grid block
825  double perf_length; // the length of the well perforation
826 
827  switch (direction) {
828  case Opm::WellCompletion::DirectionEnum::X:
829  re = std::sqrt(cubical[1] * cubical[2] / M_PI);
830  perf_length = cubical[0];
831  break;
832  case Opm::WellCompletion::DirectionEnum::Y:
833  re = std::sqrt(cubical[0] * cubical[2] / M_PI);
834  perf_length = cubical[1];
835  break;
836  case Opm::WellCompletion::DirectionEnum::Z:
837  re = std::sqrt(cubical[0] * cubical[1] / M_PI);
838  perf_length = cubical[2];
839  break;
840  default:
841  OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
842  }
843 
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);
848  }
849  }
850  }
851  }
852 
853 }
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 &param)
Constructor.
Definition: WellInterface_impl.hpp:28
WellType wellType() const
Well type, INJECTOR or PRODUCER.
Definition: WellInterface_impl.hpp:150