My Project
WellInterface_impl.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4  Copyright 2018 IRIS
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/simulators/wells/GroupState.hpp>
25 #include <opm/simulators/wells/TargetCalculator.hpp>
26 
27 #include <dune/common/version.hh>
28 
29 namespace Opm
30 {
31 
32 
33  template<typename TypeTag>
35  WellInterface(const Well& well,
36  const ParallelWellInfo& pw_info,
37  const int time_step,
38  const ModelParameters& param,
39  const RateConverterType& rate_converter,
40  const int pvtRegionIdx,
41  const int num_components,
42  const int num_phases,
43  const int index_of_well,
44  const std::vector<PerforationData>& perf_data)
45  : WellInterfaceIndices<FluidSystem,Indices,Scalar>(well,
46  pw_info,
47  time_step,
48  rate_converter,
49  pvtRegionIdx,
50  num_components,
51  num_phases,
52  index_of_well,
53  perf_data)
54  , param_(param)
55  {
56  connectionRates_.resize(this->number_of_perforations_);
57 
58  if constexpr (has_solvent || has_zFraction) {
59  if (well.isInjector()) {
60  auto injectorType = this->well_ecl_.injectorType();
61  if (injectorType == InjectorType::GAS) {
62  this->wsolvent_ = this->well_ecl_.getSolventFraction();
63  }
64  }
65  }
66  }
67 
68 
69  template<typename TypeTag>
70  void
72  init(const PhaseUsage* phase_usage_arg,
73  const std::vector<double>& /* depth_arg */,
74  const double gravity_arg,
75  const int /* num_cells */,
76  const std::vector< Scalar >& B_avg)
77  {
78  this->phase_usage_ = phase_usage_arg;
79  this->gravity_ = gravity_arg;
80  B_avg_ = B_avg;
81  }
82 
83 
84 
85 
86  template<typename TypeTag>
87  double
88  WellInterface<TypeTag>::
89  wpolymer() const
90  {
91  if constexpr (has_polymer) {
92  auto injectorType = this->well_ecl_.injectorType();
93 
94  if (injectorType == InjectorType::WATER) {
95  WellPolymerProperties polymer = this->well_ecl_.getPolymerProperties();
96  const double polymer_injection_concentration = polymer.m_polymerConcentration;
97  return polymer_injection_concentration;
98  } else {
99  // Not a water injection well => no polymer.
100  return 0.0;
101  }
102  }
103 
104  return 0.0;
105  }
106 
107 
108 
109 
110 
111  template<typename TypeTag>
112  double
113  WellInterface<TypeTag>::
114  wfoam() const
115  {
116  if constexpr (has_foam) {
117  auto injectorType = this->well_ecl_.injectorType();
118 
119  if (injectorType == InjectorType::GAS) {
120  WellFoamProperties fprop = this->well_ecl_.getFoamProperties();
121  return fprop.m_foamConcentration;
122  } else {
123  // Not a gas injection well => no foam.
124  return 0.0;
125  }
126  }
127 
128  return 0.0;
129  }
130 
131 
132 
133  template<typename TypeTag>
134  double
135  WellInterface<TypeTag>::
136  wsalt() const
137  {
138  if constexpr (has_brine) {
139  auto injectorType = this->well_ecl_.injectorType();
140 
141  if (injectorType == InjectorType::WATER) {
142  WellBrineProperties fprop = this->well_ecl_.getBrineProperties();
143  return fprop.m_saltConcentration;
144  } else {
145  // Not a water injection well => no salt (?).
146  return 0.0;
147  }
148  }
149 
150  return 0.0;
151  }
152 
153  template<typename TypeTag>
154  double
155  WellInterface<TypeTag>::
156  wmicrobes() const
157  {
158  if constexpr (has_micp) {
159  auto injectorType = this->well_ecl_.injectorType();
160 
161  if (injectorType == InjectorType::WATER) {
162  WellMICPProperties microbes = this->well_ecl_.getMICPProperties();
163  const double microbial_injection_concentration = microbes.m_microbialConcentration;
164  return microbial_injection_concentration;
165  } else {
166  // Not a water injection well => no microbes.
167  return 0.0;
168  }
169  }
170 
171  return 0.0;
172  }
173 
174  template<typename TypeTag>
175  double
176  WellInterface<TypeTag>::
177  woxygen() const
178  {
179  if constexpr (has_micp) {
180  auto injectorType = this->well_ecl_.injectorType();
181 
182  if (injectorType == InjectorType::WATER) {
183  WellMICPProperties oxygen = this->well_ecl_.getMICPProperties();
184  const double oxygen_injection_concentration = oxygen.m_oxygenConcentration;
185  return oxygen_injection_concentration;
186  } else {
187  // Not a water injection well => no oxygen.
188  return 0.0;
189  }
190  }
191 
192  return 0.0;
193  }
194 
195  // The urea injection concentration is scaled down by a factor of 10, since its value
196  // can be much bigger than 1 (not doing this slows the simulations). The
197  // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing
198  // the reactions and also when writing the output files (vtk and eclipse format, i.e.,
199  // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively).
200 
201  template<typename TypeTag>
202  double
203  WellInterface<TypeTag>::
204  wurea() const
205  {
206  if constexpr (has_micp) {
207  auto injectorType = this->well_ecl_.injectorType();
208 
209  if (injectorType == InjectorType::WATER) {
210  WellMICPProperties urea = this->well_ecl_.getMICPProperties();
211  const double urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
212  return urea_injection_concentration;
213  } else {
214  // Not a water injection well => no urea.
215  return 0.0;
216  }
217  }
218 
219  return 0.0;
220  }
221 
222  template<typename TypeTag>
223  bool
224  WellInterface<TypeTag>::
225  updateWellControl(const Simulator& ebos_simulator,
226  const IndividualOrGroup iog,
227  WellState& well_state,
228  const GroupState& group_state,
229  DeferredLogger& deferred_logger) /* const */
230  {
231  if (this->wellIsStopped()) {
232  return false;
233  }
234 
235  const auto& summaryState = ebos_simulator.vanguard().summaryState();
236  const auto& schedule = ebos_simulator.vanguard().schedule();
237  const auto& well = this->well_ecl_;
238  auto& ws = well_state.well(this->index_of_well_);
239  std::string from;
240  if (well.isInjector()) {
241  from = Well::InjectorCMode2String(ws.injection_cmode);
242  } else {
243  from = Well::ProducerCMode2String(ws.production_cmode);
244  }
245  bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
246 
247  if (oscillating) {
248  // only output frist time
249  bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
250  if (output) {
251  std::ostringstream ss;
252  ss << " The control model for well " << this->name()
253  << " is oscillating\n"
254  << " We don't allow for more than "
255  << param_.max_number_of_well_switches_
256  << " switches. The control is kept at " << from;
257  deferred_logger.info(ss.str());
258  // add one more to avoid outputting the same info again
259  this->well_control_log_.push_back(from);
260  }
261  return false;
262  }
263  bool changed = false;
264  if (iog == IndividualOrGroup::Individual) {
265  changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
266  } else if (iog == IndividualOrGroup::Group) {
267  changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
268  } else {
269  assert(iog == IndividualOrGroup::Both);
270  changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
271  }
272  Parallel::Communication cc = ebos_simulator.vanguard().grid().comm();
273  // checking whether control changed
274  if (changed) {
275  std::string to;
276  if (well.isInjector()) {
277  to = Well::InjectorCMode2String(ws.injection_cmode);
278  } else {
279  to = Well::ProducerCMode2String(ws.production_cmode);
280  }
281  std::ostringstream ss;
282  ss << " Switching control mode for well " << this->name()
283  << " from " << from
284  << " to " << to;
285  if (cc.size() > 1) {
286  ss << " on rank " << cc.rank();
287  }
288  deferred_logger.debug(ss.str());
289 
290  this->well_control_log_.push_back(from);
291  updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
292  updatePrimaryVariables(well_state, deferred_logger);
293  }
294 
295  return changed;
296  }
297 
298 
299 
300  template<typename TypeTag>
301  void
302  WellInterface<TypeTag>::
303  wellTesting(const Simulator& simulator,
304  const double simulation_time,
305  /* const */ WellState& well_state,
306  const GroupState& group_state,
307  WellTestState& well_test_state,
308  DeferredLogger& deferred_logger)
309  {
310  deferred_logger.info(" well " + this->name() + " is being tested");
311 
312  WellState well_state_copy = well_state;
313  auto& ws = well_state_copy.well(this->indexOfWell());
314 
315  updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
316  calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
317  updatePrimaryVariables(well_state_copy, deferred_logger);
318  initPrimaryVariablesEvaluation();
319 
320  WellTestState welltest_state_temp;
321 
322  bool testWell = true;
323  // if a well is closed because all completions are closed, we need to check each completion
324  // individually. We first open all completions, then we close one by one by calling updateWellTestState
325  // untill the number of closed completions do not increase anymore.
326  while (testWell) {
327  const size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
328  bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
329  if (!converged) {
330  const auto msg = fmt::format("WTEST: Well {} is not solvable (physical)", this->name());
331  deferred_logger.debug(msg);
332  return;
333  }
334 
335  updateWellOperability(simulator, well_state_copy, deferred_logger);
336  if ( !this->isOperableAndSolvable() ) {
337  const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
338  deferred_logger.debug(msg);
339  return;
340  }
341 
342  std::vector<double> potentials;
343  try {
344  computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
345  } catch (const std::exception& e) {
346  const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during testing for re-opening: ") + e.what();
347  deferred_logger.info(msg);
348  return;
349  }
350  const int np = well_state_copy.numPhases();
351  for (int p = 0; p < np; ++p) {
352  ws.well_potentials[p] = std::max(0.0, potentials[p]);
353  }
354  this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
355  this->closeCompletions(welltest_state_temp);
356 
357  // Stop testing if the well is closed or shut due to all completions shut
358  // Also check if number of completions has increased. If the number of closed completions do not increased
359  // we stop the testing.
360  // TODO: it can be tricky here, if the well is shut/closed due to other reasons
361  if ( welltest_state_temp.num_closed_wells() > 0 ||
362  (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
363  testWell = false; // this terminates the while loop
364  }
365  }
366 
367  // update wellTestState if the well test succeeds
368  if (!welltest_state_temp.well_is_closed(this->name())) {
369  well_test_state.open_well(this->name());
370 
371  std::string msg = std::string("well ") + this->name() + std::string(" is re-opened");
372  deferred_logger.info(msg);
373 
374  // also reopen completions
375  for (auto& completion : this->well_ecl_.getCompletions()) {
376  if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
377  well_test_state.open_completion(this->name(), completion.first);
378  }
379  // set the status of the well_state to open
380  ws.open();
381  well_state = well_state_copy;
382  }
383  }
384 
385 
386 
387 
388  template<typename TypeTag>
389  bool
390  WellInterface<TypeTag>::
391  iterateWellEquations(const Simulator& ebosSimulator,
392  const double dt,
393  WellState& well_state,
394  const GroupState& group_state,
395  DeferredLogger& deferred_logger)
396  {
397  const auto& summary_state = ebosSimulator.vanguard().summaryState();
398  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
399  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
400  bool converged = false;
401  try {
402  converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
403  } catch (NumericalIssue& e ) {
404  const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
405  deferred_logger.warning("INNER_ITERATION_FAILED", msg);
406  converged = false;
407  }
408  return converged;
409  }
410 
411 
412  template<typename TypeTag>
413  bool
414  WellInterface<TypeTag>::
415  solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
416  DeferredLogger& deferred_logger)
417  {
418  // keep a copy of the original well state
419  const WellState well_state0 = well_state;
420  const double dt = ebosSimulator.timeStepSize();
421  const auto& summary_state = ebosSimulator.vanguard().summaryState();
422  const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
423  if (has_thp_limit)
424  well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
425  else
426  well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
427 
428  const bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
429  if (converged) {
430  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " converged");
431  return true;
432  }
433  const int max_iter = param_.max_welleq_iter_;
434  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " failed converging in "
435  + std::to_string(max_iter) + " iterations");
436  well_state = well_state0;
437  return false;
438  }
439 
440 
441  template<typename TypeTag>
442  void
443  WellInterface<TypeTag>::
444  solveWellEquation(const Simulator& ebosSimulator,
445  WellState& well_state,
446  const GroupState& group_state,
447  DeferredLogger& deferred_logger)
448  {
449  if (!this->isOperableAndSolvable() && !this->wellIsStopped())
450  return;
451 
452  // keep a copy of the original well state
453  const WellState well_state0 = well_state;
454  const double dt = ebosSimulator.timeStepSize();
455  const bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
456  if (!converged) {
457  const int max_iter = param_.max_welleq_iter_;
458  deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
459  + std::to_string(max_iter) + " iterations");
460  well_state = well_state0;
461  }
462  }
463 
464 
465 
466  template <typename TypeTag>
467  void
468  WellInterface<TypeTag>::
469  assembleWellEq(const Simulator& ebosSimulator,
470  const double dt,
471  WellState& well_state,
472  const GroupState& group_state,
473  DeferredLogger& deferred_logger)
474  {
475  const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
476 
477  if (param_.check_well_operability_iter_)
478  checkWellOperability(ebosSimulator, well_state, deferred_logger);
479 
480  // only use inner well iterations for the first newton iterations.
481  const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
482  if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
483  this->operability_status_.solvable = true;
484  bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
485 
486  // unsolvable wells are treated as not operable and will not be solved for in this iteration.
487  if (!converged) {
488  if (param_.shut_unsolvable_wells_)
489  this->operability_status_.solvable = false;
490  }
491  }
492  if (this->operability_status_.has_negative_potentials) {
493  auto well_state_copy = well_state;
494  std::vector<double> potentials;
495  try {
496  computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
497  } catch (const std::exception& e) {
498  const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
499  deferred_logger.info(msg);
500  this->operability_status_.has_negative_potentials = true;
501  }
502  auto& ws = well_state.well(this->indexOfWell());
503  const int np = well_state.numPhases();
504  for (int p = 0; p < np; ++p) {
505  ws.well_potentials[p] = std::max(0.0, potentials[p]);
506  }
507  }
508  this->changed_to_open_this_step_ = false;
509  const bool well_operable = this->operability_status_.isOperableAndSolvable();
510 
511  if (!well_operable && old_well_operable) {
512  if (this->well_ecl_.getAutomaticShutIn()) {
513  deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
514  } else {
515  if (!this->wellIsStopped()) {
516  deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
517  this->stopWell();
518  changed_to_stopped_this_step_ = true;
519  }
520  }
521  } else if (well_operable && !old_well_operable) {
522  deferred_logger.info(" well " + this->name() + " gets REVIVED during iteration ");
523  this->openWell();
524  changed_to_stopped_this_step_ = false;
525  this->changed_to_open_this_step_ = true;
526  }
527 
528  const auto& summary_state = ebosSimulator.vanguard().summaryState();
529  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
530  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
531  assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
532  }
533 
534 
535 
536  template<typename TypeTag>
537  void
538  WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const
539  {
540  if(!this->isOperableAndSolvable() && !this->wellIsStopped())
541  return;
542 
543  for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
544  if (this->cells()[perfIdx] == cellIdx) {
545  for (int i = 0; i < RateVector::dimension; ++i) {
546  rates[i] += connectionRates_[perfIdx][i];
547  }
548  }
549  }
550  }
551 
552  template<typename TypeTag>
553  typename WellInterface<TypeTag>::Scalar
554  WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const {
555  for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
556  if (this->cells()[perfIdx] == cellIdx) {
557  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
558  return connectionRates_[perfIdx][activeCompIdx].value();
559  }
560  }
561  // this is not thread safe
562  OPM_THROW(std::invalid_argument, "The well with name " + this->name()
563  + " does not perforate cell " + std::to_string(cellIdx));
564  return 0.0;
565  }
566 
567 
568 
569 
570  template<typename TypeTag>
571  void
572  WellInterface<TypeTag>::
573  checkWellOperability(const Simulator& ebos_simulator,
574  const WellState& well_state,
575  DeferredLogger& deferred_logger)
576  {
577 
578  if (!param_.check_well_operability_) {
579  return;
580  }
581 
582  if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
583  return;
584  }
585 
586  updateWellOperability(ebos_simulator, well_state, deferred_logger);
587  }
588 
589 
590 
591  template<typename TypeTag>
592  void
593  WellInterface<TypeTag>::
594  updateWellOperability(const Simulator& ebos_simulator,
595  const WellState& well_state,
596  DeferredLogger& deferred_logger)
597  {
598  this->operability_status_.resetOperability();
599 
600  bool thp_controled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
601  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
602  bool bhp_controled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
603  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
604 
605  // Operability checking is not free
606  // Only check wells under BHP and THP control
607  bool check_thp = thp_controled || this->operability_status_.thp_limit_violated_but_not_switched;
608  if (check_thp || bhp_controled) {
609  updateIPR(ebos_simulator, deferred_logger);
610  checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
611  }
612  // we do some extra checking for wells under THP control.
613  if (thp_controled) {
614  checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
615  }
616  }
617 
618 
619  template<typename TypeTag>
620  void
621  WellInterface<TypeTag>::
622  updateWellStateWithTarget(const Simulator& ebos_simulator,
623  const GroupState& group_state,
624  WellState& well_state,
625  DeferredLogger& deferred_logger) const
626  {
627 
628  // only bhp and wellRates are used to initilize the primaryvariables for standard wells
629  const auto& well = this->well_ecl_;
630  const int well_index = this->index_of_well_;
631  auto& ws = well_state.well(well_index);
632  const auto& pu = this->phaseUsage();
633  const int np = well_state.numPhases();
634  const auto& summaryState = ebos_simulator.vanguard().summaryState();
635  const auto& schedule = ebos_simulator.vanguard().schedule();
636 
637  if (this->wellIsStopped()) {
638  for (int p = 0; p<np; ++p) {
639  ws.surface_rates[p] = 0;
640  }
641  ws.thp = 0;
642  return;
643  }
644 
645  if (this->isInjector() )
646  {
647  const auto& controls = well.injectionControls(summaryState);
648 
649  InjectorType injectorType = controls.injector_type;
650  int phasePos;
651  switch (injectorType) {
652  case InjectorType::WATER:
653  {
654  phasePos = pu.phase_pos[BlackoilPhases::Aqua];
655  break;
656  }
657  case InjectorType::OIL:
658  {
659  phasePos = pu.phase_pos[BlackoilPhases::Liquid];
660  break;
661  }
662  case InjectorType::GAS:
663  {
664  phasePos = pu.phase_pos[BlackoilPhases::Vapour];
665  break;
666  }
667  default:
668  OPM_DEFLOG_THROW(std::runtime_error, "Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
669  }
670 
671  const auto current = ws.injection_cmode;
672 
673  switch(current) {
674  case Well::InjectorCMode::RATE:
675  {
676  ws.surface_rates[phasePos] = controls.surface_rate;
677  break;
678  }
679 
680  case Well::InjectorCMode::RESV:
681  {
682  std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
683  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
684  const double coeff = convert_coeff[phasePos];
685  ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
686  break;
687  }
688 
689  case Well::InjectorCMode::THP:
690  {
691  auto rates = ws.surface_rates;
692  double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
693  ws.bhp = bhp;
694 
695  // if the total rates are negative or zero
696  // we try to provide a better intial well rate
697  // using the well potentials
698  double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
699  if (total_rate <= 0.0)
700  ws.surface_rates = ws.well_potentials;
701 
702  break;
703  }
704  case Well::InjectorCMode::BHP:
705  {
706  ws.bhp = controls.bhp_limit;
707  double total_rate = 0.0;
708  for (int p = 0; p<np; ++p) {
709  total_rate += ws.surface_rates[p];
710  }
711  // if the total rates are negative or zero
712  // we try to provide a better intial well rate
713  // using the well potentials
714  if (total_rate <= 0.0)
715  ws.surface_rates = ws.well_potentials;
716 
717  break;
718  }
719  case Well::InjectorCMode::GRUP:
720  {
721  assert(well.isAvailableForGroupControl());
722  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
723  const double efficiencyFactor = well.getEfficiencyFactor();
724  std::optional<double> target =
725  this->getGroupInjectionTargetRate(group,
726  well_state,
727  group_state,
728  schedule,
729  summaryState,
730  injectorType,
731  efficiencyFactor,
732  deferred_logger);
733  if (target)
734  ws.surface_rates[phasePos] = *target;
735  break;
736  }
737  case Well::InjectorCMode::CMODE_UNDEFINED:
738  {
739  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger );
740  }
741 
742  }
743  }
744  //Producer
745  else
746  {
747  const auto current = ws.production_cmode;
748  const auto& controls = well.productionControls(summaryState);
749  switch (current) {
750  case Well::ProducerCMode::ORAT:
751  {
752  double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
753  // for trivial rates or opposite direction we don't just scale the rates
754  // but use either the potentials or the mobility ratio to initial the well rates
755  if (current_rate > 0.0) {
756  for (int p = 0; p<np; ++p) {
757  ws.surface_rates[p] *= controls.oil_rate/current_rate;
758  }
759  } else {
760  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
761  double control_fraction = fractions[pu.phase_pos[Oil]];
762  if (control_fraction != 0.0) {
763  for (int p = 0; p<np; ++p) {
764  ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
765  }
766  }
767  }
768  break;
769  }
770  case Well::ProducerCMode::WRAT:
771  {
772  double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
773  // for trivial rates or opposite direction we don't just scale the rates
774  // but use either the potentials or the mobility ratio to initial the well rates
775  if (current_rate > 0.0) {
776  for (int p = 0; p<np; ++p) {
777  ws.surface_rates[p] *= controls.water_rate/current_rate;
778  }
779  } else {
780  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
781  double control_fraction = fractions[pu.phase_pos[Water]];
782  if (control_fraction != 0.0) {
783  for (int p = 0; p<np; ++p) {
784  ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
785  }
786  }
787  }
788  break;
789  }
790  case Well::ProducerCMode::GRAT:
791  {
792  double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
793  // or trivial rates or opposite direction we don't just scale the rates
794  // but use either the potentials or the mobility ratio to initial the well rates
795  if (current_rate > 0.0) {
796  for (int p = 0; p<np; ++p) {
797  ws.surface_rates[p] *= controls.gas_rate/current_rate;
798  }
799  } else {
800  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
801  double control_fraction = fractions[pu.phase_pos[Gas]];
802  if (control_fraction != 0.0) {
803  for (int p = 0; p<np; ++p) {
804  ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
805  }
806  }
807  }
808 
809  break;
810 
811  }
812  case Well::ProducerCMode::LRAT:
813  {
814  double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
815  - ws.surface_rates[ pu.phase_pos[Oil] ];
816  // or trivial rates or opposite direction we don't just scale the rates
817  // but use either the potentials or the mobility ratio to initial the well rates
818  if (current_rate > 0.0) {
819  for (int p = 0; p<np; ++p) {
820  ws.surface_rates[p] *= controls.liquid_rate/current_rate;
821  }
822  } else {
823  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
824  double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
825  if (control_fraction != 0.0) {
826  for (int p = 0; p<np; ++p) {
827  ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
828  }
829  }
830  }
831  break;
832  }
833  case Well::ProducerCMode::CRAT:
834  {
835  OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << this->name(), deferred_logger);
836  }
837  case Well::ProducerCMode::RESV:
838  {
839  std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
840  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
841  double total_res_rate = 0.0;
842  for (int p = 0; p<np; ++p) {
843  total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
844  }
845  if (controls.prediction_mode) {
846  // or trivial rates or opposite direction we don't just scale the rates
847  // but use either the potentials or the mobility ratio to initial the well rates
848  if (total_res_rate > 0.0) {
849  for (int p = 0; p<np; ++p) {
850  ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
851  }
852  } else {
853  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
854  for (int p = 0; p<np; ++p) {
855  ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
856  }
857  }
858  } else {
859  std::vector<double> hrates(this->number_of_phases_,0.);
860  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
861  hrates[pu.phase_pos[Water]] = controls.water_rate;
862  }
863  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
864  hrates[pu.phase_pos[Oil]] = controls.oil_rate;
865  }
866  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
867  hrates[pu.phase_pos[Gas]] = controls.gas_rate;
868  }
869  std::vector<double> hrates_resv(this->number_of_phases_,0.);
870  this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
871  double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
872  // or trivial rates or opposite direction we don't just scale the rates
873  // but use either the potentials or the mobility ratio to initial the well rates
874  if (total_res_rate > 0.0) {
875  for (int p = 0; p<np; ++p) {
876  ws.surface_rates[p] *= target/total_res_rate;
877  }
878  } else {
879  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
880  for (int p = 0; p<np; ++p) {
881  ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
882  }
883  }
884 
885  }
886  break;
887  }
888  case Well::ProducerCMode::BHP:
889  {
890  ws.bhp = controls.bhp_limit;
891  double total_rate = 0.0;
892  for (int p = 0; p<np; ++p) {
893  total_rate -= ws.surface_rates[p];
894  }
895  // if the total rates are negative or zero
896  // we try to provide a better intial well rate
897  // using the well potentials
898  if (total_rate <= 0.0){
899  for (int p = 0; p<np; ++p) {
900  ws.surface_rates[p] = -ws.well_potentials[p];
901  }
902  }
903  break;
904  }
905  case Well::ProducerCMode::THP:
906  {
907  auto rates = ws.surface_rates;
908  this->adaptRatesForVFP(rates);
909  double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
910  ws.bhp = bhp;
911 
912  // if the total rates are negative or zero
913  // we try to provide a better intial well rate
914  // using the well potentials
915  double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
916  if (total_rate <= 0.0){
917  for (int p = 0; p<np; ++p) {
918  ws.surface_rates[p] = -ws.well_potentials[p];
919  }
920  }
921  break;
922  }
923  case Well::ProducerCMode::GRUP:
924  {
925  assert(well.isAvailableForGroupControl());
926  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
927  const double efficiencyFactor = well.getEfficiencyFactor();
928  double scale = this->getGroupProductionTargetRate(group,
929  well_state,
930  group_state,
931  schedule,
932  summaryState,
933  efficiencyFactor);
934 
935  // we don't want to scale with zero and get zero rates.
936  if (scale > 0) {
937  for (int p = 0; p<np; ++p) {
938  ws.surface_rates[p] *= scale;
939  }
940  }
941  break;
942  }
943  case Well::ProducerCMode::CMODE_UNDEFINED:
944  case Well::ProducerCMode::NONE:
945  {
946  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name() , deferred_logger);
947  }
948 
949  break;
950  } // end of switch
951  }
952  }
953 
954  template<typename TypeTag>
955  std::vector<double>
956  WellInterface<TypeTag>::
957  initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const
958  {
959  const int np = this->number_of_phases_;
960  std::vector<double> scaling_factor(np);
961  const auto& ws = well_state.well(this->index_of_well_);
962 
963  double total_potentials = 0.0;
964  for (int p = 0; p<np; ++p) {
965  total_potentials += ws.well_potentials[p];
966  }
967  if (total_potentials > 0) {
968  for (int p = 0; p<np; ++p) {
969  scaling_factor[p] = ws.well_potentials[p] / total_potentials;
970  }
971  return scaling_factor;
972  }
973  // if we don't have any potentials we weight it using the mobilites
974  // We only need approximation so we don't bother with the vapporized oil and dissolved gas
975  double total_tw = 0;
976  const int nperf = this->number_of_perforations_;
977  for (int perf = 0; perf < nperf; ++perf) {
978  total_tw += this->well_index_[perf];
979  }
980  for (int perf = 0; perf < nperf; ++perf) {
981  const int cell_idx = this->well_cells_[perf];
982  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
983  const auto& fs = intQuants.fluidState();
984  const double well_tw_fraction = this->well_index_[perf] / total_tw;
985  double total_mobility = 0.0;
986  for (int p = 0; p < np; ++p) {
987  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
988  total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
989  }
990  for (int p = 0; p < np; ++p) {
991  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
992  scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
993  }
994  }
995  return scaling_factor;
996  }
997 
998 
999 
1000  template <typename TypeTag>
1001  void
1003  updateWellStateRates(const Simulator& ebosSimulator,
1004  WellState& well_state,
1005  DeferredLogger& deferred_logger) const
1006  {
1007  // Check if the rates of this well only are single-phase, do nothing
1008  // if more than one nonzero rate.
1009  auto& ws = well_state.well(this->index_of_well_);
1010  int nonzero_rate_index = -1;
1011  for (int p = 0; p < this->number_of_phases_; ++p) {
1012  if (ws.surface_rates[p] != 0.0) {
1013  if (nonzero_rate_index == -1) {
1014  nonzero_rate_index = p;
1015  } else {
1016  // More than one nonzero rate.
1017  return;
1018  }
1019  }
1020  }
1021 
1022  // Calculate the rates that follow from the current primary variables.
1023  std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
1024 
1025  if (nonzero_rate_index == -1) {
1026  // No nonzero rates.
1027  // Use the computed rate directly
1028  for (int p = 0; p < this->number_of_phases_; ++p) {
1029  ws.surface_rates[p] = well_q_s[this->flowPhaseToEbosCompIdx(p)];
1030  }
1031  return;
1032  }
1033 
1034  // Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
1035  const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1036  const int comp_idx_nz = this->flowPhaseToEbosCompIdx(nonzero_rate_index);
1037  for (int p = 0; p < this->number_of_phases_; ++p) {
1038  if (p != nonzero_rate_index) {
1039  const int comp_idx = this->flowPhaseToEbosCompIdx(p);
1040  double& rate = ws.surface_rates[p];
1041  rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
1042  }
1043  }
1044  }
1045  template<typename TypeTag>
1046  typename WellInterface<TypeTag>::Eval
1047  WellInterface<TypeTag>::getPerfCellPressure(const typename WellInterface<TypeTag>::FluidState& fs) const
1048  {
1049  Eval pressure;
1050  if (Indices::oilEnabled) {
1051  pressure = fs.pressure(FluidSystem::oilPhaseIdx);
1052  } else {
1053  if (Indices::waterEnabled) {
1054  pressure = fs.pressure(FluidSystem::waterPhaseIdx);
1055  } else {
1056  pressure = fs.pressure(FluidSystem::gasPhaseIdx);
1057  }
1058  }
1059  return pressure;
1060  }
1061 } // namespace Opm
Definition: DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:252
Definition: WellInterfaceIndices.hpp:35
Definition: WellInterface.hpp:72
void updateWellStateRates(const Simulator &ebosSimulator, WellState &well_state, DeferredLogger &deferred_logger) const
Modify the well_state's rates if there is only one nonzero rate.
Definition: WellInterface_impl.hpp:1003
WellInterface(const Well &well, const ParallelWellInfo &pw_info, const int time_step, const ModelParameters &param, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_components, const int num_phases, const int index_of_well, const std::vector< PerforationData > &perf_data)
Constructor.
Definition: WellInterface_impl.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParametersEbos.hpp:327
Definition: BlackoilPhases.hpp:46