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  const bool changed_to_open_this_step)
78  {
79  this->phase_usage_ = phase_usage_arg;
80  this->gravity_ = gravity_arg;
81  B_avg_ = B_avg;
82  this->changed_to_open_this_step_ = changed_to_open_this_step;
83  }
84 
85 
86 
87 
88  template<typename TypeTag>
89  double
90  WellInterface<TypeTag>::
91  wpolymer() const
92  {
93  if constexpr (has_polymer) {
94  auto injectorType = this->well_ecl_.injectorType();
95 
96  if (injectorType == InjectorType::WATER) {
97  WellPolymerProperties polymer = this->well_ecl_.getPolymerProperties();
98  const double polymer_injection_concentration = polymer.m_polymerConcentration;
99  return polymer_injection_concentration;
100  } else {
101  // Not a water injection well => no polymer.
102  return 0.0;
103  }
104  }
105 
106  return 0.0;
107  }
108 
109 
110 
111 
112 
113  template<typename TypeTag>
114  double
115  WellInterface<TypeTag>::
116  wfoam() const
117  {
118  if constexpr (has_foam) {
119  auto injectorType = this->well_ecl_.injectorType();
120 
121  if (injectorType == InjectorType::GAS) {
122  WellFoamProperties fprop = this->well_ecl_.getFoamProperties();
123  return fprop.m_foamConcentration;
124  } else {
125  // Not a gas injection well => no foam.
126  return 0.0;
127  }
128  }
129 
130  return 0.0;
131  }
132 
133 
134 
135  template<typename TypeTag>
136  double
137  WellInterface<TypeTag>::
138  wsalt() const
139  {
140  if constexpr (has_brine) {
141  auto injectorType = this->well_ecl_.injectorType();
142 
143  if (injectorType == InjectorType::WATER) {
144  WellBrineProperties fprop = this->well_ecl_.getBrineProperties();
145  return fprop.m_saltConcentration;
146  } else {
147  // Not a water injection well => no salt (?).
148  return 0.0;
149  }
150  }
151 
152  return 0.0;
153  }
154 
155  template<typename TypeTag>
156  double
157  WellInterface<TypeTag>::
158  wmicrobes() const
159  {
160  if constexpr (has_micp) {
161  auto injectorType = this->well_ecl_.injectorType();
162 
163  if (injectorType == InjectorType::WATER) {
164  WellMICPProperties microbes = this->well_ecl_.getMICPProperties();
165  const double microbial_injection_concentration = microbes.m_microbialConcentration;
166  return microbial_injection_concentration;
167  } else {
168  // Not a water injection well => no microbes.
169  return 0.0;
170  }
171  }
172 
173  return 0.0;
174  }
175 
176  template<typename TypeTag>
177  double
178  WellInterface<TypeTag>::
179  woxygen() const
180  {
181  if constexpr (has_micp) {
182  auto injectorType = this->well_ecl_.injectorType();
183 
184  if (injectorType == InjectorType::WATER) {
185  WellMICPProperties oxygen = this->well_ecl_.getMICPProperties();
186  const double oxygen_injection_concentration = oxygen.m_oxygenConcentration;
187  return oxygen_injection_concentration;
188  } else {
189  // Not a water injection well => no oxygen.
190  return 0.0;
191  }
192  }
193 
194  return 0.0;
195  }
196 
197  // The urea injection concentration is scaled down by a factor of 10, since its value
198  // can be much bigger than 1 (not doing this slows the simulations). The
199  // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing
200  // the reactions and also when writing the output files (vtk and eclipse format, i.e.,
201  // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively).
202 
203  template<typename TypeTag>
204  double
205  WellInterface<TypeTag>::
206  wurea() const
207  {
208  if constexpr (has_micp) {
209  auto injectorType = this->well_ecl_.injectorType();
210 
211  if (injectorType == InjectorType::WATER) {
212  WellMICPProperties urea = this->well_ecl_.getMICPProperties();
213  const double urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
214  return urea_injection_concentration;
215  } else {
216  // Not a water injection well => no urea.
217  return 0.0;
218  }
219  }
220 
221  return 0.0;
222  }
223 
224  template<typename TypeTag>
225  bool
226  WellInterface<TypeTag>::
227  updateWellControl(const Simulator& ebos_simulator,
228  const IndividualOrGroup iog,
229  WellState& well_state,
230  const GroupState& group_state,
231  DeferredLogger& deferred_logger) /* const */
232  {
233  if (this->wellIsStopped()) {
234  return false;
235  }
236 
237  const auto& summaryState = ebos_simulator.vanguard().summaryState();
238  const auto& schedule = ebos_simulator.vanguard().schedule();
239  const auto& well = this->well_ecl_;
240  auto& ws = well_state.well(this->index_of_well_);
241  std::string from;
242  if (well.isInjector()) {
243  from = Well::InjectorCMode2String(ws.injection_cmode);
244  } else {
245  from = Well::ProducerCMode2String(ws.production_cmode);
246  }
247  bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
248 
249  if (oscillating) {
250  // only output frist time
251  bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
252  if (output) {
253  std::ostringstream ss;
254  ss << " The control model for well " << this->name()
255  << " is oscillating\n"
256  << " We don't allow for more than "
257  << param_.max_number_of_well_switches_
258  << " switches. The control is kept at " << from;
259  deferred_logger.info(ss.str());
260  // add one more to avoid outputting the same info again
261  this->well_control_log_.push_back(from);
262  }
263  return false;
264  }
265  bool changed = false;
266  if (iog == IndividualOrGroup::Individual) {
267  changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
268  } else if (iog == IndividualOrGroup::Group) {
269  changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
270  } else {
271  assert(iog == IndividualOrGroup::Both);
272  changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
273  }
274  Parallel::Communication cc = ebos_simulator.vanguard().grid().comm();
275  // checking whether control changed
276  if (changed) {
277  std::string to;
278  if (well.isInjector()) {
279  to = Well::InjectorCMode2String(ws.injection_cmode);
280  } else {
281  to = Well::ProducerCMode2String(ws.production_cmode);
282  }
283  std::ostringstream ss;
284  ss << " Switching control mode for well " << this->name()
285  << " from " << from
286  << " to " << to;
287  if (cc.size() > 1) {
288  ss << " on rank " << cc.rank();
289  }
290  deferred_logger.debug(ss.str());
291 
292  this->well_control_log_.push_back(from);
293  updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
294  updatePrimaryVariables(well_state, deferred_logger);
295  }
296 
297  return changed;
298  }
299 
300 
301 
302  template<typename TypeTag>
303  void
304  WellInterface<TypeTag>::
305  wellTesting(const Simulator& simulator,
306  const double simulation_time,
307  /* const */ WellState& well_state,
308  const GroupState& group_state,
309  WellTestState& well_test_state,
310  DeferredLogger& deferred_logger)
311  {
312  deferred_logger.info(" well " + this->name() + " is being tested");
313 
314  WellState well_state_copy = well_state;
315  auto& ws = well_state_copy.well(this->indexOfWell());
316 
317  updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
318  calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
319  updatePrimaryVariables(well_state_copy, deferred_logger);
320  initPrimaryVariablesEvaluation();
321 
322  if (this->isProducer()) {
323  gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
324  }
325 
326  WellTestState welltest_state_temp;
327 
328  bool testWell = true;
329  // if a well is closed because all completions are closed, we need to check each completion
330  // individually. We first open all completions, then we close one by one by calling updateWellTestState
331  // untill the number of closed completions do not increase anymore.
332  while (testWell) {
333  const size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
334  bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
335  if (!converged) {
336  const auto msg = fmt::format("WTEST: Well {} is not solvable (physical)", this->name());
337  deferred_logger.debug(msg);
338  return;
339  }
340 
341 
342  updateWellOperability(simulator, well_state_copy, deferred_logger);
343  if ( !this->isOperableAndSolvable() ) {
344  const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
345  deferred_logger.debug(msg);
346  return;
347  }
348 
349  std::vector<double> potentials;
350  try {
351  computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
352  } catch (const std::exception& e) {
353  const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during testing for re-opening: ") + e.what();
354  deferred_logger.info(msg);
355  return;
356  }
357  const int np = well_state_copy.numPhases();
358  for (int p = 0; p < np; ++p) {
359  ws.well_potentials[p] = std::max(0.0, potentials[p]);
360  }
361  this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
362  this->closeCompletions(welltest_state_temp);
363 
364  // Stop testing if the well is closed or shut due to all completions shut
365  // Also check if number of completions has increased. If the number of closed completions do not increased
366  // we stop the testing.
367  // TODO: it can be tricky here, if the well is shut/closed due to other reasons
368  if ( welltest_state_temp.num_closed_wells() > 0 ||
369  (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
370  testWell = false; // this terminates the while loop
371  }
372  }
373 
374  // update wellTestState if the well test succeeds
375  if (!welltest_state_temp.well_is_closed(this->name())) {
376  well_test_state.open_well(this->name());
377 
378  std::string msg = std::string("well ") + this->name() + std::string(" is re-opened");
379  deferred_logger.info(msg);
380 
381  // also reopen completions
382  for (auto& completion : this->well_ecl_.getCompletions()) {
383  if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
384  well_test_state.open_completion(this->name(), completion.first);
385  }
386  // set the status of the well_state to open
387  ws.open();
388  well_state = well_state_copy;
389  }
390  }
391 
392 
393 
394 
395  template<typename TypeTag>
396  bool
397  WellInterface<TypeTag>::
398  iterateWellEquations(const Simulator& ebosSimulator,
399  const double dt,
400  WellState& well_state,
401  const GroupState& group_state,
402  DeferredLogger& deferred_logger)
403  {
404  const auto& summary_state = ebosSimulator.vanguard().summaryState();
405  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
406  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
407  bool converged = false;
408  try {
409  converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
410  } catch (NumericalIssue& e ) {
411  const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
412  deferred_logger.warning("INNER_ITERATION_FAILED", msg);
413  converged = false;
414  }
415  return converged;
416  }
417 
418 
419  template<typename TypeTag>
420  bool
421  WellInterface<TypeTag>::
422  solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
423  DeferredLogger& deferred_logger)
424  {
425  // keep a copy of the original well state
426  const WellState well_state0 = well_state;
427  const double dt = ebosSimulator.timeStepSize();
428  const auto& summary_state = ebosSimulator.vanguard().summaryState();
429  const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
430  if (has_thp_limit)
431  well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
432  else
433  well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
434 
435  const bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
436  if (converged) {
437  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " converged");
438  return true;
439  }
440  const int max_iter = param_.max_welleq_iter_;
441  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " failed converging in "
442  + std::to_string(max_iter) + " iterations");
443  well_state = well_state0;
444  return false;
445  }
446 
447 
448  template<typename TypeTag>
449  void
450  WellInterface<TypeTag>::
451  solveWellEquation(const Simulator& ebosSimulator,
452  WellState& well_state,
453  const GroupState& group_state,
454  DeferredLogger& deferred_logger)
455  {
456  if (!this->isOperableAndSolvable() && !this->wellIsStopped())
457  return;
458 
459  // keep a copy of the original well state
460  const WellState well_state0 = well_state;
461  const double dt = ebosSimulator.timeStepSize();
462  bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
463 
464  // Newly opened wells with THP control sometimes struggles to
465  // converge due to bad initial guess. Or due to the simple fact
466  // that the well needs to change to another control.
467  // We therefore try to solve the well with BHP control to get
468  // an better initial guess.
469  // If the well is supposed to operate under THP control
470  // "updateWellControl" will switch it back to THP later.
471  if (!converged) {
472  auto& ws = well_state.well(this->indexOfWell());
473  bool thp_control = false;
474  if (this->well_ecl_.isInjector()) {
475  thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
476  if (thp_control) {
477  ws.injection_cmode = Well::InjectorCMode::BHP;
478  this->well_control_log_.push_back(Well::InjectorCMode2String(Well::InjectorCMode::THP));
479  }
480  } else {
481  thp_control = ws.production_cmode == Well::ProducerCMode::THP;
482  if (thp_control) {
483  ws.production_cmode = Well::ProducerCMode::BHP;
484  this->well_control_log_.push_back(Well::ProducerCMode2String(Well::ProducerCMode::THP));
485  }
486  }
487  if (thp_control) {
488  const std::string msg = std::string("The newly opened well ") + this->name()
489  + std::string(" with THP control did not converge during inner iterations, we try again with bhp control");
490  deferred_logger.debug(msg);
491  converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
492  }
493  }
494 
495  if (!converged) {
496  const int max_iter = param_.max_welleq_iter_;
497  deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
498  + std::to_string(max_iter) + " iterations");
499  well_state = well_state0;
500  }
501  }
502 
503 
504 
505  template <typename TypeTag>
506  void
507  WellInterface<TypeTag>::
508  assembleWellEq(const Simulator& ebosSimulator,
509  const double dt,
510  WellState& well_state,
511  const GroupState& group_state,
512  DeferredLogger& deferred_logger)
513  {
514  const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
515 
516  if (param_.check_well_operability_iter_)
517  checkWellOperability(ebosSimulator, well_state, deferred_logger);
518 
519  // only use inner well iterations for the first newton iterations.
520  const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
521  if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
522  this->operability_status_.solvable = true;
523  bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
524 
525  // unsolvable wells are treated as not operable and will not be solved for in this iteration.
526  if (!converged) {
527  if (param_.shut_unsolvable_wells_)
528  this->operability_status_.solvable = false;
529  }
530  }
531  if (this->operability_status_.has_negative_potentials) {
532  auto well_state_copy = well_state;
533  std::vector<double> potentials;
534  try {
535  computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
536  } catch (const std::exception& e) {
537  const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
538  deferred_logger.info(msg);
539  this->operability_status_.has_negative_potentials = true;
540  }
541  auto& ws = well_state.well(this->indexOfWell());
542  const int np = well_state.numPhases();
543  for (int p = 0; p < np; ++p) {
544  ws.well_potentials[p] = std::max(0.0, potentials[p]);
545  }
546  }
547  this->changed_to_open_this_step_ = false;
548  const bool well_operable = this->operability_status_.isOperableAndSolvable();
549 
550  if (!well_operable && old_well_operable) {
551  if (this->well_ecl_.getAutomaticShutIn()) {
552  deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
553  } else {
554  if (!this->wellIsStopped()) {
555  deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
556  this->stopWell();
557  changed_to_stopped_this_step_ = true;
558  }
559  }
560  } else if (well_operable && !old_well_operable) {
561  deferred_logger.info(" well " + this->name() + " gets REVIVED during iteration ");
562  this->openWell();
563  changed_to_stopped_this_step_ = false;
564  this->changed_to_open_this_step_ = true;
565  }
566 
567  const auto& summary_state = ebosSimulator.vanguard().summaryState();
568  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
569  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
570  assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
571  }
572 
573  template<typename TypeTag>
574  bool
575  WellInterface<TypeTag>::isPressureControlled(const WellState& well_state) const
576  {
577  bool thp_controlled_well = false;
578  bool bhp_controlled_well = false;
579  const auto& ws = well_state.well(this->index_of_well_);
580  if (this->isInjector()) {
581  const Well::InjectorCMode& current = ws.injection_cmode;
582  if (current == Well::InjectorCMode::THP) {
583  thp_controlled_well = true;
584  }
585  if (current == Well::InjectorCMode::BHP) {
586  bhp_controlled_well = true;
587  }
588  } else {
589  const Well::ProducerCMode& current = ws.production_cmode;
590  if (current == Well::ProducerCMode::THP) {
591  thp_controlled_well = true;
592  }
593  if (current == Well::ProducerCMode::BHP) {
594  bhp_controlled_well = true;
595  }
596  }
597  bool ispressureControlled = (bhp_controlled_well || thp_controlled_well);
598  return ispressureControlled;
599  }
600 
601  template<typename TypeTag>
602  void
603  WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const
604  {
605  if(!this->isOperableAndSolvable() && !this->wellIsStopped())
606  return;
607 
608  for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
609  if (this->cells()[perfIdx] == cellIdx) {
610  for (int i = 0; i < RateVector::dimension; ++i) {
611  rates[i] += connectionRates_[perfIdx][i];
612  }
613  }
614  }
615  }
616 
617  template<typename TypeTag>
618  typename WellInterface<TypeTag>::Scalar
619  WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const {
620  for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
621  if (this->cells()[perfIdx] == cellIdx) {
622  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
623  return connectionRates_[perfIdx][activeCompIdx].value();
624  }
625  }
626  // this is not thread safe
627  OPM_THROW(std::invalid_argument, "The well with name " + this->name()
628  + " does not perforate cell " + std::to_string(cellIdx));
629  return 0.0;
630  }
631 
632 
633 
634 
635  template<typename TypeTag>
636  void
637  WellInterface<TypeTag>::
638  checkWellOperability(const Simulator& ebos_simulator,
639  const WellState& well_state,
640  DeferredLogger& deferred_logger)
641  {
642 
643  if (!param_.check_well_operability_) {
644  return;
645  }
646 
647  if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
648  return;
649  }
650 
651  updateWellOperability(ebos_simulator, well_state, deferred_logger);
652  if (!this->operability_status_.isOperableAndSolvable()) {
653  this->operability_status_.use_vfpexplicit = true;
654  deferred_logger.debug("EXPLICIT_LOOKUP_VFP",
655  "well not operable, trying with explicit vfp lookup: " + this->name());
656  updateWellOperability(ebos_simulator, well_state, deferred_logger);
657  }
658  }
659 
660  template<typename TypeTag>
661  void
662  WellInterface<TypeTag>::
663  gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
664  WellState& well_state,
665  DeferredLogger& deferred_logger)
666  {
667  const auto& summary_state = ebos_simulator.vanguard().summaryState();
668  const auto& well_name = this->name();
669  if (!this->wellHasTHPConstraints(summary_state)) {
670  const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
671  deferred_logger.info(msg);
672  return;
673  }
674  const auto& well_ecl = this->wellEcl();
675  const auto& schedule = ebos_simulator.vanguard().schedule();
676  auto report_step_idx = ebos_simulator.episodeIndex();
677  const auto& glo = schedule.glo(report_step_idx);
678  if (!glo.has_well(well_name)) {
679  const std::string msg = fmt::format(
680  "GLIFT WTEST: Well {} : Gas Lift not activated: "
681  "WLIFTOPT is probably missing. Skipping.", well_name);
682  deferred_logger.info(msg);
683  return;
684  }
685  const auto& gl_well = glo.well(well_name);
686  auto& max_alq_optional = gl_well.max_rate();
687  double max_alq;
688  if (max_alq_optional) {
689  max_alq = *max_alq_optional;
690  }
691  else {
692  const auto& controls = well_ecl.productionControls(summary_state);
693  const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
694  const auto& alq_values = table.getALQAxis();
695  max_alq = alq_values.back();
696  }
697  well_state.setALQ(well_name, max_alq);
698  const std::string msg = fmt::format(
699  "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
700  well_name, max_alq);
701  deferred_logger.info(msg);
702  }
703 
704  template<typename TypeTag>
705  void
706  WellInterface<TypeTag>::
707  updateWellOperability(const Simulator& ebos_simulator,
708  const WellState& well_state,
709  DeferredLogger& deferred_logger)
710  {
711  this->operability_status_.resetOperability();
712 
713  bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
714  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
715  bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
716  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
717 
718  // Operability checking is not free
719  // Only check wells under BHP and THP control
720  bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
721  if (check_thp || bhp_controlled) {
722  updateIPR(ebos_simulator, deferred_logger);
723  checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
724  }
725  // we do some extra checking for wells under THP control.
726  if (check_thp) {
727  checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
728  }
729  }
730 
731 
732  template<typename TypeTag>
733  void
734  WellInterface<TypeTag>::
735  updateWellStateWithTarget(const Simulator& ebos_simulator,
736  const GroupState& group_state,
737  WellState& well_state,
738  DeferredLogger& deferred_logger) const
739  {
740 
741  // only bhp and wellRates are used to initilize the primaryvariables for standard wells
742  const auto& well = this->well_ecl_;
743  const int well_index = this->index_of_well_;
744  auto& ws = well_state.well(well_index);
745  const auto& pu = this->phaseUsage();
746  const int np = well_state.numPhases();
747  const auto& summaryState = ebos_simulator.vanguard().summaryState();
748  const auto& schedule = ebos_simulator.vanguard().schedule();
749 
750  if (this->wellIsStopped()) {
751  for (int p = 0; p<np; ++p) {
752  ws.surface_rates[p] = 0;
753  }
754  ws.thp = 0;
755  return;
756  }
757 
758  if (this->isInjector() )
759  {
760  const auto& controls = well.injectionControls(summaryState);
761 
762  InjectorType injectorType = controls.injector_type;
763  int phasePos;
764  switch (injectorType) {
765  case InjectorType::WATER:
766  {
767  phasePos = pu.phase_pos[BlackoilPhases::Aqua];
768  break;
769  }
770  case InjectorType::OIL:
771  {
772  phasePos = pu.phase_pos[BlackoilPhases::Liquid];
773  break;
774  }
775  case InjectorType::GAS:
776  {
777  phasePos = pu.phase_pos[BlackoilPhases::Vapour];
778  break;
779  }
780  default:
781  OPM_DEFLOG_THROW(std::runtime_error, "Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
782  }
783 
784  const auto current = ws.injection_cmode;
785 
786  switch(current) {
787  case Well::InjectorCMode::RATE:
788  {
789  ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
790  if(this->rsRvInj() > 0) {
791  if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
792  ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
793  } else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
794  ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
795  } else {
796  OPM_DEFLOG_THROW(std::runtime_error, "Expected OIL or GAS as type for injectors when RS/RV (item 10) is non-zero " + this->name(), deferred_logger );
797  }
798  }
799  break;
800  }
801 
802  case Well::InjectorCMode::RESV:
803  {
804  std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
805  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
806  const double coeff = convert_coeff[phasePos];
807  ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
808  break;
809  }
810 
811  case Well::InjectorCMode::THP:
812  {
813  auto rates = ws.surface_rates;
814  double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
815  ws.bhp = bhp;
816  ws.thp = this->getTHPConstraint(summaryState);
817 
818  // if the total rates are negative or zero
819  // we try to provide a better intial well rate
820  // using the well potentials
821  double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
822  if (total_rate <= 0.0)
823  ws.surface_rates = ws.well_potentials;
824 
825  break;
826  }
827  case Well::InjectorCMode::BHP:
828  {
829  ws.bhp = controls.bhp_limit;
830  double total_rate = 0.0;
831  for (int p = 0; p<np; ++p) {
832  total_rate += ws.surface_rates[p];
833  }
834  // if the total rates are negative or zero
835  // we try to provide a better intial well rate
836  // using the well potentials
837  if (total_rate <= 0.0)
838  ws.surface_rates = ws.well_potentials;
839 
840  break;
841  }
842  case Well::InjectorCMode::GRUP:
843  {
844  assert(well.isAvailableForGroupControl());
845  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
846  const double efficiencyFactor = well.getEfficiencyFactor();
847  std::optional<double> target =
848  this->getGroupInjectionTargetRate(group,
849  well_state,
850  group_state,
851  schedule,
852  summaryState,
853  injectorType,
854  efficiencyFactor,
855  deferred_logger);
856  if (target)
857  ws.surface_rates[phasePos] = *target;
858  break;
859  }
860  case Well::InjectorCMode::CMODE_UNDEFINED:
861  {
862  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger );
863  }
864 
865  }
866  // for wells with zero injection rate, if we assign exactly zero rate,
867  // we will have to assume some trivial composition in the wellbore.
868  // here, we use some small value (about 0.01 m^3/day ~= 1.e-7) to initialize
869  // the zero rate target, then we can use to retain the composition information
870  // within the wellbore from the previous result, and hopefully it is a good
871  // initial guess for the zero rate target.
872  ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
873  }
874  //Producer
875  else
876  {
877  const auto current = ws.production_cmode;
878  const auto& controls = well.productionControls(summaryState);
879  switch (current) {
880  case Well::ProducerCMode::ORAT:
881  {
882  double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
883  // for trivial rates or opposite direction we don't just scale the rates
884  // but use either the potentials or the mobility ratio to initial the well rates
885  if (current_rate > 0.0) {
886  for (int p = 0; p<np; ++p) {
887  ws.surface_rates[p] *= controls.oil_rate/current_rate;
888  }
889  } else {
890  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
891  double control_fraction = fractions[pu.phase_pos[Oil]];
892  if (control_fraction != 0.0) {
893  for (int p = 0; p<np; ++p) {
894  ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
895  }
896  }
897  }
898  break;
899  }
900  case Well::ProducerCMode::WRAT:
901  {
902  double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
903  // for trivial rates or opposite direction we don't just scale the rates
904  // but use either the potentials or the mobility ratio to initial the well rates
905  if (current_rate > 0.0) {
906  for (int p = 0; p<np; ++p) {
907  ws.surface_rates[p] *= controls.water_rate/current_rate;
908  }
909  } else {
910  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
911  double control_fraction = fractions[pu.phase_pos[Water]];
912  if (control_fraction != 0.0) {
913  for (int p = 0; p<np; ++p) {
914  ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
915  }
916  }
917  }
918  break;
919  }
920  case Well::ProducerCMode::GRAT:
921  {
922  double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
923  // or trivial rates or opposite direction we don't just scale the rates
924  // but use either the potentials or the mobility ratio to initial the well rates
925  if (current_rate > 0.0) {
926  for (int p = 0; p<np; ++p) {
927  ws.surface_rates[p] *= controls.gas_rate/current_rate;
928  }
929  } else {
930  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
931  double control_fraction = fractions[pu.phase_pos[Gas]];
932  if (control_fraction != 0.0) {
933  for (int p = 0; p<np; ++p) {
934  ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
935  }
936  }
937  }
938 
939  break;
940 
941  }
942  case Well::ProducerCMode::LRAT:
943  {
944  double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
945  - ws.surface_rates[ pu.phase_pos[Oil] ];
946  // or trivial rates or opposite direction we don't just scale the rates
947  // but use either the potentials or the mobility ratio to initial the well rates
948  if (current_rate > 0.0) {
949  for (int p = 0; p<np; ++p) {
950  ws.surface_rates[p] *= controls.liquid_rate/current_rate;
951  }
952  } else {
953  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
954  double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
955  if (control_fraction != 0.0) {
956  for (int p = 0; p<np; ++p) {
957  ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
958  }
959  }
960  }
961  break;
962  }
963  case Well::ProducerCMode::CRAT:
964  {
965  OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << this->name(), deferred_logger);
966  }
967  case Well::ProducerCMode::RESV:
968  {
969  std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
970  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
971  double total_res_rate = 0.0;
972  for (int p = 0; p<np; ++p) {
973  total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
974  }
975  if (controls.prediction_mode) {
976  // or trivial rates or opposite direction we don't just scale the rates
977  // but use either the potentials or the mobility ratio to initial the well rates
978  if (total_res_rate > 0.0) {
979  for (int p = 0; p<np; ++p) {
980  ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
981  }
982  } else {
983  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
984  for (int p = 0; p<np; ++p) {
985  ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
986  }
987  }
988  } else {
989  std::vector<double> hrates(this->number_of_phases_,0.);
990  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
991  hrates[pu.phase_pos[Water]] = controls.water_rate;
992  }
993  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
994  hrates[pu.phase_pos[Oil]] = controls.oil_rate;
995  }
996  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
997  hrates[pu.phase_pos[Gas]] = controls.gas_rate;
998  }
999  std::vector<double> hrates_resv(this->number_of_phases_,0.);
1000  this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
1001  double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1002  // or trivial rates or opposite direction we don't just scale the rates
1003  // but use either the potentials or the mobility ratio to initial the well rates
1004  if (total_res_rate > 0.0) {
1005  for (int p = 0; p<np; ++p) {
1006  ws.surface_rates[p] *= target/total_res_rate;
1007  }
1008  } else {
1009  const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1010  for (int p = 0; p<np; ++p) {
1011  ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1012  }
1013  }
1014 
1015  }
1016  break;
1017  }
1018  case Well::ProducerCMode::BHP:
1019  {
1020  ws.bhp = controls.bhp_limit;
1021  double total_rate = 0.0;
1022  for (int p = 0; p<np; ++p) {
1023  total_rate -= ws.surface_rates[p];
1024  }
1025  // if the total rates are negative or zero
1026  // we try to provide a better intial well rate
1027  // using the well potentials
1028  if (total_rate <= 0.0){
1029  for (int p = 0; p<np; ++p) {
1030  ws.surface_rates[p] = -ws.well_potentials[p];
1031  }
1032  }
1033  break;
1034  }
1035  case Well::ProducerCMode::THP:
1036  {
1037  auto rates = ws.surface_rates;
1038  this->adaptRatesForVFP(rates);
1039  double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1040  ws.bhp = bhp;
1041  ws.thp = this->getTHPConstraint(summaryState);
1042 
1043  // if the total rates are negative or zero
1044  // we try to provide a better intial well rate
1045  // using the well potentials
1046  double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1047  if (total_rate <= 0.0){
1048  for (int p = 0; p<np; ++p) {
1049  ws.surface_rates[p] = -ws.well_potentials[p];
1050  }
1051  }
1052  break;
1053  }
1054  case Well::ProducerCMode::GRUP:
1055  {
1056  assert(well.isAvailableForGroupControl());
1057  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1058  const double efficiencyFactor = well.getEfficiencyFactor();
1059  double scale = this->getGroupProductionTargetRate(group,
1060  well_state,
1061  group_state,
1062  schedule,
1063  summaryState,
1064  efficiencyFactor);
1065 
1066  // we don't want to scale with zero and get zero rates.
1067  if (scale > 0) {
1068  for (int p = 0; p<np; ++p) {
1069  ws.surface_rates[p] *= scale;
1070  }
1071  ws.trivial_target = false;
1072  } else {
1073  ws.trivial_target = true;
1074  }
1075  break;
1076  }
1077  case Well::ProducerCMode::CMODE_UNDEFINED:
1078  case Well::ProducerCMode::NONE:
1079  {
1080  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name() , deferred_logger);
1081  }
1082 
1083  break;
1084  } // end of switch
1085  }
1086  }
1087 
1088  template<typename TypeTag>
1089  std::vector<double>
1090  WellInterface<TypeTag>::
1091  initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const
1092  {
1093  const int np = this->number_of_phases_;
1094  std::vector<double> scaling_factor(np);
1095  const auto& ws = well_state.well(this->index_of_well_);
1096 
1097  double total_potentials = 0.0;
1098  for (int p = 0; p<np; ++p) {
1099  total_potentials += ws.well_potentials[p];
1100  }
1101  if (total_potentials > 0) {
1102  for (int p = 0; p<np; ++p) {
1103  scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1104  }
1105  return scaling_factor;
1106  }
1107  // if we don't have any potentials we weight it using the mobilites
1108  // We only need approximation so we don't bother with the vapporized oil and dissolved gas
1109  double total_tw = 0;
1110  const int nperf = this->number_of_perforations_;
1111  for (int perf = 0; perf < nperf; ++perf) {
1112  total_tw += this->well_index_[perf];
1113  }
1114  for (int perf = 0; perf < nperf; ++perf) {
1115  const int cell_idx = this->well_cells_[perf];
1116  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1117  const auto& fs = intQuants.fluidState();
1118  const double well_tw_fraction = this->well_index_[perf] / total_tw;
1119  double total_mobility = 0.0;
1120  for (int p = 0; p < np; ++p) {
1121  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1122  total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1123  }
1124  for (int p = 0; p < np; ++p) {
1125  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1126  scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1127  }
1128  }
1129  return scaling_factor;
1130  }
1131 
1132 
1133 
1134  template <typename TypeTag>
1135  void
1137  updateWellStateRates(const Simulator& ebosSimulator,
1138  WellState& well_state,
1139  DeferredLogger& deferred_logger) const
1140  {
1141  // Check if the rates of this well only are single-phase, do nothing
1142  // if more than one nonzero rate.
1143  auto& ws = well_state.well(this->index_of_well_);
1144  int nonzero_rate_index = -1;
1145  const double floating_point_error_epsilon = 1e-14;
1146  for (int p = 0; p < this->number_of_phases_; ++p) {
1147  if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1148  if (nonzero_rate_index == -1) {
1149  nonzero_rate_index = p;
1150  } else {
1151  // More than one nonzero rate.
1152  return;
1153  }
1154  }
1155  }
1156 
1157  // Calculate the rates that follow from the current primary variables.
1158  std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
1159 
1160  if (nonzero_rate_index == -1) {
1161  // No nonzero rates.
1162  // Use the computed rate directly
1163  for (int p = 0; p < this->number_of_phases_; ++p) {
1164  ws.surface_rates[p] = well_q_s[this->flowPhaseToEbosCompIdx(p)];
1165  }
1166  return;
1167  }
1168 
1169  // Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
1170  const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1171  const int comp_idx_nz = this->flowPhaseToEbosCompIdx(nonzero_rate_index);
1172  for (int p = 0; p < this->number_of_phases_; ++p) {
1173  if (p != nonzero_rate_index) {
1174  const int comp_idx = this->flowPhaseToEbosCompIdx(p);
1175  double& rate = ws.surface_rates[p];
1176  rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
1177  }
1178  }
1179  }
1180  template<typename TypeTag>
1181  typename WellInterface<TypeTag>::Eval
1182  WellInterface<TypeTag>::getPerfCellPressure(const typename WellInterface<TypeTag>::FluidState& fs) const
1183  {
1184  Eval pressure;
1185  if (Indices::oilEnabled) {
1186  pressure = fs.pressure(FluidSystem::oilPhaseIdx);
1187  } else {
1188  if (Indices::waterEnabled) {
1189  pressure = fs.pressure(FluidSystem::waterPhaseIdx);
1190  } else {
1191  pressure = fs.pressure(FluidSystem::gasPhaseIdx);
1192  }
1193  }
1194  return pressure;
1195  }
1196 } // namespace Opm
Definition: DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:243
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:1137
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:37
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParametersEbos.hpp:327
Definition: BlackoilPhases.hpp:46