opm-simulators
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 #ifndef OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
23 #define OPM_WELLINTERFACE_IMPL_HEADER_INCLUDED
24 
25 // Improve IDE experience
26 #ifndef OPM_WELLINTERFACE_HEADER_INCLUDED
27 #include <config.h>
28 #include <opm/simulators/wells/WellInterface.hpp>
29 #endif
30 
31 #include <opm/common/Exceptions.hpp>
32 
33 #include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
34 #include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
35 
36 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
37 
38 #include <opm/simulators/wells/GroupState.hpp>
39 #include <opm/simulators/wells/TargetCalculator.hpp>
40 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
41 #include <opm/simulators/wells/WellHelpers.hpp>
42 
43 #include <dune/common/version.hh>
44 
45 #include <algorithm>
46 #include <cassert>
47 #include <cstddef>
48 #include <numbers>
49 #include <utility>
50 
51 #include <fmt/format.h>
52 
53 namespace Opm
54 {
55 
56 
57  template<typename TypeTag>
58  WellInterface<TypeTag>::
59  WellInterface(const Well& well,
60  const ParallelWellInfo<Scalar>& pw_info,
61  const int time_step,
62  const ModelParameters& param,
63  const RateConverterType& rate_converter,
64  const int pvtRegionIdx,
65  const int num_conservation_quantities,
66  const int num_phases,
67  const int index_of_well,
68  const std::vector<PerforationData<Scalar>>& perf_data)
69  : WellInterfaceIndices<FluidSystem,Indices>(well,
70  pw_info,
71  time_step,
72  param,
73  rate_converter,
74  pvtRegionIdx,
75  num_conservation_quantities,
76  num_phases,
77  index_of_well,
78  perf_data)
79  {
80  connectionRates_.resize(this->number_of_local_perforations_);
81 
82  if constexpr (has_solvent || has_zFraction) {
83  if (well.isInjector()) {
84  auto injectorType = this->well_ecl_.injectorType();
85  if (injectorType == InjectorType::GAS) {
86  this->wsolvent_ = this->well_ecl_.getSolventFraction();
87  }
88  }
89  }
90  }
91 
92 
93  template<typename TypeTag>
94  void
95  WellInterface<TypeTag>::
96  init(const std::vector<Scalar>& /* depth_arg */,
97  const Scalar gravity_arg,
98  const std::vector<Scalar>& B_avg,
99  const bool changed_to_open_this_step)
100  {
101  this->gravity_ = gravity_arg;
102  B_avg_ = B_avg;
103  this->changed_to_open_this_step_ = changed_to_open_this_step;
104  }
105 
106 
107 
108 
109  template<typename TypeTag>
110  typename WellInterface<TypeTag>::Scalar
111  WellInterface<TypeTag>::
112  wpolymer() const
113  {
114  if constexpr (has_polymer) {
115  return this->wpolymer_();
116  }
117 
118  return 0.0;
119  }
120 
121 
122 
123 
124 
125  template<typename TypeTag>
126  typename WellInterface<TypeTag>::Scalar
127  WellInterface<TypeTag>::
128  wfoam() const
129  {
130  if constexpr (has_foam) {
131  return this->wfoam_();
132  }
133 
134  return 0.0;
135  }
136 
137 
138 
139  template<typename TypeTag>
140  typename WellInterface<TypeTag>::Scalar
141  WellInterface<TypeTag>::
142  wsalt() const
143  {
144  if constexpr (has_brine) {
145  return this->wsalt_();
146  }
147 
148  return 0.0;
149  }
150 
151  template<typename TypeTag>
152  typename WellInterface<TypeTag>::Scalar
153  WellInterface<TypeTag>::
154  wmicrobes() const
155  {
156  if constexpr (has_micp) {
157  return this->wmicrobes_();
158  }
159 
160  return 0.0;
161  }
162 
163  template<typename TypeTag>
164  typename WellInterface<TypeTag>::Scalar
165  WellInterface<TypeTag>::
166  woxygen() const
167  {
168  if constexpr (has_micp) {
169  return this->woxygen_();
170  }
171 
172  return 0.0;
173  }
174 
175  template<typename TypeTag>
176  typename WellInterface<TypeTag>::Scalar
177  WellInterface<TypeTag>::
178  wurea() const
179  {
180  if constexpr (has_micp) {
181  return this->wurea_();
182  }
183 
184  return 0.0;
185  }
186 
187  template<typename TypeTag>
188  bool
189  WellInterface<TypeTag>::
190  updateWellControl(const Simulator& simulator,
191  const IndividualOrGroup iog,
192  const GroupStateHelperType& groupStateHelper,
193  WellStateType& well_state) /* const */
194  {
195  auto& deferred_logger = groupStateHelper.deferredLogger();
196  OPM_TIMEFUNCTION();
197  if (stoppedOrZeroRateTarget(groupStateHelper)) {
198  return false;
199  }
200 
201  const auto& summaryState = simulator.vanguard().summaryState();
202  const auto& schedule = simulator.vanguard().schedule();
203  const auto& well = this->well_ecl_;
204  auto& ws = well_state.well(this->index_of_well_);
205  std::string from;
206  bool is_grup = false;
207  if (well.isInjector()) {
208  from = WellInjectorCMode2String(ws.injection_cmode);
209  is_grup = ws.injection_cmode == Well::InjectorCMode::GRUP;
210  } else {
211  from = WellProducerCMode2String(ws.production_cmode);
212  is_grup = ws.production_cmode == Well::ProducerCMode::GRUP;
213  }
214 
215  const int episodeIdx = simulator.episodeIndex();
216  const auto& iterCtx = simulator.problem().iterationContext();
217  const int nupcol = schedule[episodeIdx].nupcol();
218  const bool oscillating =
219  std::ranges::count(this->well_control_log_, from) >= this->param_.max_number_of_well_switches_;
220  if (oscillating && !is_grup) { // we would like to avoid ending up as GRUP
221  // only output first time
222  const bool output =
223  std::ranges::count(this->well_control_log_, from) == this->param_.max_number_of_well_switches_;
224  if (output) {
225  const auto msg = fmt::format(" The control mode for well {} is oscillating. \n"
226  "We don't allow for more than {} switches after NUPCOL iterations. (NUPCOL = {}) \n"
227  "The control is kept at {}.",
228  this->name(), this->param_.max_number_of_well_switches_, nupcol, from);
229  deferred_logger.info(msg);
230  // add one more to avoid outputting the same info again
231  this->well_control_log_.push_back(from);
232  }
233  return false;
234  }
235  bool changed = false;
236  if (iog == IndividualOrGroup::Individual) {
237  changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
238  } else if (iog == IndividualOrGroup::Group) {
239  changed = this->checkGroupConstraints(
240  groupStateHelper, schedule, summaryState, true, well_state
241  );
242  } else {
243  assert(iog == IndividualOrGroup::Both);
244  changed = this->checkConstraints(groupStateHelper, schedule, summaryState, well_state);
245  }
246  Parallel::Communication cc = simulator.vanguard().grid().comm();
247  // checking whether control changed
248  if (changed) {
249  std::string to;
250  if (well.isInjector()) {
251  to = WellInjectorCMode2String(ws.injection_cmode);
252  } else {
253  to = WellProducerCMode2String(ws.production_cmode);
254  }
255  std::ostringstream ss;
256  ss << " Switching control mode for well " << this->name()
257  << " from " << from
258  << " to " << to;
259  if (iterCtx.inLocalSolve()) {
260  ss << " (NLDD domain solve)";
261  }
262  if (cc.size() > 1) {
263  ss << " on rank " << cc.rank();
264  }
265  deferred_logger.debug(ss.str());
266 
267  // We always store the current control as it is used for output
268  // and only after iteration >= nupcol
269  // we log all switches to check if the well controls oscillates.
270  // Skip logging during NLDD local solves to avoid exhausting the
271  // global oscillation budget with provisional domain-level switches.
272  if (!iterCtx.inLocalSolve()) {
273  if (!iterCtx.withinNupcol(nupcol) || this->well_control_log_.empty()) {
274  this->well_control_log_.push_back(from);
275  }
276  }
277  updateWellStateWithTarget(simulator, groupStateHelper, well_state);
278  updatePrimaryVariables(groupStateHelper);
279  }
280 
281  return changed;
282  }
283 
284  template<typename TypeTag>
285  bool
286  WellInterface<TypeTag>::
287  updateWellControlAndStatusLocalIteration(const Simulator& simulator,
288  const GroupStateHelperType& groupStateHelper,
289  const Well::InjectionControls& inj_controls,
290  const Well::ProductionControls& prod_controls,
291  const Scalar wqTotal,
292  WellStateType& well_state,
293  const bool fixed_control,
294  const bool fixed_status,
295  const bool solving_with_zero_rate)
296  {
297  OPM_TIMEFUNCTION();
298  auto& deferred_logger = groupStateHelper.deferredLogger();
299  const auto& summary_state = simulator.vanguard().summaryState();
300  const auto& schedule = simulator.vanguard().schedule();
301  auto& ws = well_state.well(this->index_of_well_);
302  std::string from;
303  if (this->isInjector()) {
304  from = WellInjectorCMode2String(ws.injection_cmode);
305  } else {
306  from = WellProducerCMode2String(ws.production_cmode);
307  }
308  const bool oscillating =
309  std::ranges::count(this->well_control_log_, from) >= this->param_.max_number_of_well_switches_;
310 
311  if (oscillating || this->wellUnderZeroRateTarget(groupStateHelper) || !(well_state.well(this->index_of_well_).status == WellStatus::OPEN)) {
312  return false;
313  }
314 
315  const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
316  if (!this->wellIsStopped()){
317  if (wqTotal*sgn <= 0.0 && !fixed_status){
318  this->stopWell();
319  return true;
320  } else {
321  bool changed = false;
322  if (!fixed_control) {
323  // When solving_with_zero_rate=true, fixed_control=true, so this block should never
324  // be entered.
325  if (solving_with_zero_rate) {
326  OPM_DEFLOG_THROW(std::runtime_error, fmt::format(
327  "Well {}: solving_with_zero_rate should not be true when fixed_control is false",
328  this->name()), deferred_logger);
329  }
330 
331  // Changing to group controls here may lead to inconsistencies in the group handling which in turn
332  // may result in excessive back and forth switching. However, we currently allow this by default.
333  // The switch check_group_constraints_inner_well_iterations_ is a temporary solution.
334  const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
335  prod_controls.hasControl(Well::ProducerCMode::GRUP);
336  bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
337  if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
338  changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
339  }
340  if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
341  changed = changed || this->checkGroupConstraints(
342  groupStateHelper, schedule, summary_state, false, well_state
343  );
344  }
345 
346  if (changed) {
347  const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
348  ws.production_cmode == Well::ProducerCMode::THP;
349  if (thp_controlled){
350  ws.thp = this->getTHPConstraint(summary_state);
351  } else {
352  // don't call for thp since this might trigger additional local solve
353  updateWellStateWithTarget(simulator, groupStateHelper, well_state);
354  }
355  updatePrimaryVariables(groupStateHelper);
356  }
357  }
358  return changed;
359  }
360  } else if (!fixed_status){
361  // well is stopped, check if current bhp allows reopening
362  const Scalar bhp = well_state.well(this->index_of_well_).bhp;
363  Scalar prod_limit = prod_controls.bhp_limit;
364  Scalar inj_limit = inj_controls.bhp_limit;
365  const bool has_thp = this->wellHasTHPConstraints(summary_state);
366  if (has_thp){
367  std::vector<Scalar> rates(this->num_conservation_quantities_);
368  if (this->isInjector()){
369  const Scalar bhp_thp = WellBhpThpCalculator(*this).
370  calculateBhpFromThp(well_state, rates,
371  this->well_ecl_,
372  summary_state,
373  this->getRefDensity(),
374  deferred_logger);
375  inj_limit = std::min(bhp_thp, static_cast<Scalar>(inj_controls.bhp_limit));
376  } else {
377  // if the well can operate, it must at least be able to produce
378  // at the lowest bhp of the bhp-curve (explicit fractions)
379  const Scalar bhp_min = WellBhpThpCalculator(*this).
380  calculateMinimumBhpFromThp(well_state,
381  this->well_ecl_,
382  summary_state,
383  this->getRefDensity());
384  prod_limit = std::max(bhp_min, static_cast<Scalar>(prod_controls.bhp_limit));
385  }
386  }
387  const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
388  if (bhp_diff > 0){
389  this->openWell();
390  well_state.well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
391  if (has_thp) {
392  well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
393  }
394  return true;
395  } else {
396  return false;
397  }
398  } else {
399  return false;
400  }
401  }
402 
403  template<typename TypeTag>
404  void
405  WellInterface<TypeTag>::
406  wellTesting(const Simulator& simulator,
407  const double simulation_time,
408  const GroupStateHelperType& groupStateHelper,
409  WellStateType& well_state,
410  WellTestState& well_test_state,
411  GLiftEclWells& ecl_well_map,
412  std::map<std::string, double>& open_times)
413  {
414  OPM_TIMEFUNCTION();
415  auto& deferred_logger = groupStateHelper.deferredLogger();
416  const auto& group_state = groupStateHelper.groupState();
417  deferred_logger.info(" well " + this->name() + " is being tested");
418 
419  GroupStateHelperType groupStateHelper_copy = groupStateHelper;
420  WellStateType well_state_copy = well_state;
421  // Ensure that groupStateHelper uses well_state_copy as WellState for the well testing
422  // and the guard ensures that the original well state is restored at scope exit, i.e. at
423  // the end of this function.
424  auto guard = groupStateHelper_copy.pushWellState(well_state_copy);
425  auto& ws = well_state_copy.well(this->indexOfWell());
426 
427  const auto& summary_state = simulator.vanguard().summaryState();
428  const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
429  if (this->isProducer()) {
430  ws.production_cmode = has_thp_limit ? Well::ProducerCMode::THP : Well::ProducerCMode::BHP;
431  } else {
432  ws.injection_cmode = has_thp_limit ? Well::InjectorCMode::THP : Well::InjectorCMode::BHP;
433  }
434  // We test the well as an open well during the well testing
435  ws.open();
436 
437  scaleSegmentRatesAndPressure(well_state_copy);
438  calculateExplicitQuantities(simulator, groupStateHelper_copy);
439  updatePrimaryVariables(groupStateHelper_copy);
440 
441  if (this->isProducer()) {
442  const auto& schedule = simulator.vanguard().schedule();
443  const auto report_step = simulator.episodeIndex();
444  const auto& glo = schedule.glo(report_step);
445  if (glo.active()) {
446  gliftBeginTimeStepWellTestUpdateALQ(simulator,
447  well_state_copy,
448  group_state,
449  ecl_well_map,
450  deferred_logger);
451  }
452  }
453 
454  WellTestState welltest_state_temp;
455 
456  bool testWell = true;
457  // if a well is closed because all completions are closed, we need to check each completion
458  // individually. We first open all completions, then we close one by one by calling updateWellTestState
459  // untill the number of closed completions do not increase anymore.
460  while (testWell) {
461  const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
462  bool converged = solveWellForTesting(simulator, groupStateHelper_copy, well_state_copy);
463  if (!converged) {
464  const auto msg = fmt::format("WTEST: Well {} is not solvable (physical)", this->name());
465  deferred_logger.debug(msg);
466  return;
467  }
468 
469 
470  updateWellOperability(simulator, well_state_copy, groupStateHelper_copy);
471  if ( !this->isOperableAndSolvable() ) {
472  const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
473  deferred_logger.debug(msg);
474  return;
475  }
476  std::vector<Scalar> potentials;
477  try {
478  computeWellPotentials(simulator, well_state_copy, groupStateHelper_copy, potentials);
479  } catch (const std::exception& e) {
480  const std::string msg = fmt::format("well {}: computeWellPotentials() "
481  "failed during testing for re-opening: ",
482  this->name(), e.what());
483  deferred_logger.info(msg);
484  return;
485  }
486  const int np = well_state_copy.numPhases();
487  for (int p = 0; p < np; ++p) {
488  ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
489  }
490  const bool under_zero_target = this->wellUnderZeroGroupRateTarget(groupStateHelper_copy);
491  this->updateWellTestState(well_state_copy.well(this->indexOfWell()),
492  simulation_time,
493  /*writeMessageToOPMLog=*/ false,
494  under_zero_target,
495  welltest_state_temp,
496  deferred_logger);
497  this->closeCompletions(welltest_state_temp);
498 
499  // Stop testing if the well is closed or shut due to all completions shut
500  // Also check if number of completions has increased. If the number of closed completions do not increased
501  // we stop the testing.
502  // TODO: it can be tricky here, if the well is shut/closed due to other reasons
503  if ( welltest_state_temp.num_closed_wells() > 0 ||
504  (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
505  testWell = false; // this terminates the while loop
506  }
507  }
508 
509  // update wellTestState if the well test succeeds
510  if (!welltest_state_temp.well_is_closed(this->name())) {
511  well_test_state.open_well(this->name());
512 
513  std::string msg = std::string("well ") + this->name() + std::string(" is re-opened");
514  deferred_logger.info(msg);
515 
516  // also reopen completions
517  for (const auto& completion : this->well_ecl_.getCompletions()) {
518  if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
519  well_test_state.open_completion(this->name(), completion.first);
520  }
521  well_state = well_state_copy;
522  open_times.try_emplace(this->name(), well_test_state.lastTestTime(this->name()));
523  }
524  }
525 
526 
527 
528 
529  template<typename TypeTag>
530  bool
531  WellInterface<TypeTag>::
532  iterateWellEquations(const Simulator& simulator,
533  const double dt,
534  const GroupStateHelperType& groupStateHelper,
535  WellStateType& well_state)
536  {
537  OPM_TIMEFUNCTION();
538  auto& deferred_logger = groupStateHelper.deferredLogger();
539 
540  const auto& summary_state = simulator.vanguard().summaryState();
541  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
542  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
543  const auto& ws = well_state.well(this->indexOfWell());
544  const auto pmode_orig = ws.production_cmode;
545  const auto imode_orig = ws.injection_cmode;
546  bool converged = false;
547  try {
548  // TODO: the following two functions will be refactored to be one to reduce the code duplication
549  if (!this->param_.local_well_solver_control_switching_){
550  converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state);
551  } else {
552  if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.well(this->index_of_well_).status == WellStatus::OPEN)) {
553  converged = solveWellWithOperabilityCheck(
554  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state
555  );
556  } else {
557  converged = this->iterateWellEqWithSwitching(
558  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state,
559  /*fixed_control=*/false, /*fixed_status=*/false, /*solving_with_zero_rate=*/false
560  );
561  }
562  }
563 
564  } catch (NumericalProblem& e ) {
565  const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
566  deferred_logger.warning("INNER_ITERATION_FAILED", msg);
567  converged = false;
568  }
569  if (converged) {
570  // Add debug info for switched controls
571  if (ws.production_cmode != pmode_orig || ws.injection_cmode != imode_orig) {
572  std::string from,to;
573  if (this->isInjector()) {
574  from = WellInjectorCMode2String(imode_orig);
575  to = WellInjectorCMode2String(ws.injection_cmode);
576  } else {
577  from = WellProducerCMode2String(pmode_orig);
578  to = WellProducerCMode2String(ws.production_cmode);
579  }
580  const auto msg = fmt::format(" Well {} switched from {} to {} during local solve", this->name(), from, to);
581  deferred_logger.debug(msg);
582  const int episodeIdx = simulator.episodeIndex();
583  const auto& iterCtx = simulator.problem().iterationContext();
584  const auto& schedule = simulator.vanguard().schedule();
585  const int nupcol = schedule[episodeIdx].nupcol();
586  // We always store the current control as it is used for output
587  // and only after iteration >= nupcol
588  // we log all switches to check if the well controls oscillates
589  if (!iterCtx.withinNupcol(nupcol) || this->well_control_log_.empty()) {
590  this->well_control_log_.push_back(from);
591  }
592  }
593  }
594  // Add debug info for problematic group targets
595  if (this->isProducer() && ws.production_cmode == Well::ProducerCMode::GRUP && ws.use_group_target_fallback) {
596  assert(ws.group_target && ws.group_target_fallback);
597  const std::string cmode = Group::ProductionCMode2String(ws.group_target->production_cmode);
598  const std::string cmode_fallback = Group::ProductionCMode2String(ws.group_target_fallback->production_cmode);
599  const auto msg = fmt::format(" Well {} was solved using group target fallback mode {} as current group mode {} was not feasible.",
600  this->name(), cmode_fallback, cmode);
601  deferred_logger.debug(msg);
602  }
603 
604  return converged;
605  }
606 
607  template<typename TypeTag>
608  bool
609  WellInterface<TypeTag>::
610  solveWellWithOperabilityCheck(const Simulator& simulator,
611  const double dt,
612  const Well::InjectionControls& inj_controls,
613  const Well::ProductionControls& prod_controls,
614  const GroupStateHelperType& groupStateHelper,
615  WellStateType& well_state)
616  {
617  OPM_TIMEFUNCTION();
618  auto& deferred_logger = groupStateHelper.deferredLogger();
619 
620  const auto& summary_state = simulator.vanguard().summaryState();
621  bool converged = true;
622  auto& ws = well_state.well(this->index_of_well_);
623  // if well is stopped, check if we can reopen with explicit fraction
624  if (this->wellIsStopped()) {
625  this->openWell();
626  const bool use_vfpexplicit = this->operability_status_.use_vfpexplicit;
627  this->operability_status_.use_vfpexplicit = true;
628  auto bhp_target = estimateOperableBhp(simulator, dt, groupStateHelper, summary_state, well_state);
629  if (!bhp_target.has_value()) {
630  // no intersection with ipr
631  const auto msg = fmt::format("estimateOperableBhp: Did not find operable BHP for well {}", this->name());
632  deferred_logger.debug(msg);
633  // well can't operate using explicit fractions stop the well
634  // solve with zero rates
635  converged = solveWellWithZeroRate(simulator, dt, groupStateHelper, well_state);
636  this->stopWell();
637  this->operability_status_.can_obtain_bhp_with_thp_limit = false;
638  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
639  return converged;
640  } else {
641  // solve well with the estimated target bhp (or limit)
642  ws.thp = this->getTHPConstraint(summary_state);
643  const Scalar bhp = std::max(bhp_target.value(),
644  static_cast<Scalar>(prod_controls.bhp_limit));
645  solveWellWithBhp(simulator, dt, bhp, groupStateHelper, well_state);
646  this->operability_status_.use_vfpexplicit = use_vfpexplicit;
647  }
648  }
649  // solve well-equation
650  converged = this->iterateWellEqWithSwitching(
651  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state,
652  /*fixed_control=*/false, /*fixed_status=*/false, /*solving_with_zero_rate=*/false
653  );
654 
655 
656  const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
657  // check stability of solution under thp-control
658  if (converged && !stoppedOrZeroRateTarget(groupStateHelper) && isThp) {
659  auto rates = well_state.well(this->index_of_well_).surface_rates;
660  this->adaptRatesForVFP(rates);
661  this->updateIPRImplicit(simulator, groupStateHelper, well_state);
662  bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state);
663  if (!is_stable) {
664  // solution converged to an unstable point!
665  this->operability_status_.use_vfpexplicit = true;
666  auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
667  // if we find an intersection with a sufficiently lower bhp, re-solve equations
668  const Scalar reltol = 1e-3;
669  const Scalar cur_bhp = ws.bhp;
670  if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
671  const auto msg = fmt::format("Well {} converged to an unstable solution, re-solving", this->name());
672  deferred_logger.debug(msg);
673  solveWellWithBhp(
674  simulator, dt, bhp_stable.value(), groupStateHelper, well_state
675  );
676  // re-solve with hopefully good initial guess
677  ws.thp = this->getTHPConstraint(summary_state);
678  converged = this->iterateWellEqWithSwitching(
679  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state,
680  /*fixed_control=*/false, /*fixed_status=*/false, /*solving_with_zero_rate=*/false
681  );
682  }
683  }
684  }
685 
686  if (!converged) {
687  // Well did not converge, switch to explicit fractions
688  this->operability_status_.use_vfpexplicit = true;
689  this->openWell();
690  auto bhp_target = estimateOperableBhp(
691  simulator, dt, groupStateHelper, summary_state, well_state
692  );
693  if (!bhp_target.has_value()) {
694  // solve with zero rate
695  // well can't operate using explicit fractions stop the well
696  converged = solveWellWithZeroRate(simulator, dt, groupStateHelper, well_state);
697  this->stopWell();
698  this->operability_status_.can_obtain_bhp_with_thp_limit = false;
699  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
700  return converged;
701  } else {
702  // solve well with the estimated target bhp (or limit)
703  const Scalar bhp = std::max(bhp_target.value(),
704  static_cast<Scalar>(prod_controls.bhp_limit));
705  solveWellWithBhp(
706  simulator, dt, bhp, groupStateHelper, well_state
707  );
708  ws.thp = this->getTHPConstraint(summary_state);
709  const auto msg = fmt::format("Well {} did not converge, re-solving with explicit fractions for VFP caculations.", this->name());
710  deferred_logger.debug(msg);
711  converged = this->iterateWellEqWithSwitching(simulator, dt,
712  inj_controls,
713  prod_controls,
714  groupStateHelper,
715  well_state,
716  /*fixed_control=*/false,
717  /*fixed_status=*/false,
718  /*solving_with_zero_rate=*/false);
719  }
720  }
721  // update operability
722  this->operability_status_.can_obtain_bhp_with_thp_limit = !this->wellIsStopped();
723  this->operability_status_.obey_thp_limit_under_bhp_limit = !this->wellIsStopped();
724  return converged;
725  }
726 
727  template<typename TypeTag>
728  std::optional<typename WellInterface<TypeTag>::Scalar>
729  WellInterface<TypeTag>::
730  estimateOperableBhp(const Simulator& simulator,
731  const double dt,
732  const GroupStateHelperType& groupStateHelper,
733  const SummaryState& summary_state,
734  WellStateType& well_state)
735  {
736  if (!this->wellHasTHPConstraints(summary_state)) {
737  const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summary_state);
738  const bool converged = solveWellWithBhp(
739  simulator, dt, bhp_limit, groupStateHelper, well_state
740  );
741  if (!converged || this->wellIsStopped()) {
742  return std::nullopt;
743  }
744 
745  return bhp_limit;
746  }
747  OPM_TIMEFUNCTION();
748  // Given an unconverged well or closed well, estimate an operable bhp (if any)
749  // Get minimal bhp from vfp-curve
750  Scalar bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
751  // Solve
752  const bool converged = solveWellWithBhp(
753  simulator, dt, bhp_min, groupStateHelper, well_state
754  );
755  if (!converged || this->wellIsStopped()) {
756  return std::nullopt;
757  }
758  this->updateIPRImplicit(simulator, groupStateHelper, well_state);
759  auto rates = well_state.well(this->index_of_well_).surface_rates;
760  this->adaptRatesForVFP(rates);
761  return WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
762  }
763 
764  template<typename TypeTag>
765  bool
766  WellInterface<TypeTag>::
767  solveWellWithBhp(const Simulator& simulator,
768  const double dt,
769  const Scalar bhp,
770  const GroupStateHelperType& groupStateHelper,
771  WellStateType& well_state)
772  {
773  OPM_TIMEFUNCTION();
774 
775  // Solve a well using single bhp-constraint (but close if not operable under this)
776  auto group_state = GroupState<Scalar>(); // empty group
777  GroupStateHelperType groupStateHelper_copy = groupStateHelper;
778  // Ensure that groupStateHelper_copy uses the empty group state as GroupState for iterateWellEqWithSwitching()
779  // and the guard ensures that the original group state is restored at scope exit, i.e. at
780  // the end of this function.
781  auto group_guard = groupStateHelper_copy.pushGroupState(group_state);
782 
783  auto inj_controls = Well::InjectionControls(0);
784  auto prod_controls = Well::ProductionControls(0);
785  auto& ws = well_state.well(this->index_of_well_);
786  auto cmode_inj = ws.injection_cmode;
787  auto cmode_prod = ws.production_cmode;
788  if (this->isInjector()) {
789  inj_controls.addControl(Well::InjectorCMode::BHP);
790  inj_controls.bhp_limit = bhp;
791  inj_controls.cmode = Well::InjectorCMode::BHP;
792  ws.injection_cmode = Well::InjectorCMode::BHP;
793  } else {
794  prod_controls.addControl(Well::ProducerCMode::BHP);
795  prod_controls.bhp_limit = bhp;
796  prod_controls.cmode = Well::ProducerCMode::BHP;
797  ws.production_cmode = Well::ProducerCMode::BHP;
798  }
799  // update well-state
800  ws.bhp = bhp;
801  // solve
802  const bool converged = this->iterateWellEqWithSwitching(
803  simulator, dt, inj_controls, prod_controls, groupStateHelper_copy,
804  well_state,
805  /*fixed_control=*/true,
806  /*fixed_status=*/false,
807  /*solving_with_zero_rate=*/false
808  );
809  ws.injection_cmode = cmode_inj;
810  ws.production_cmode = cmode_prod;
811  return converged;
812  }
813 
814  template<typename TypeTag>
815  bool
816  WellInterface<TypeTag>::
817  solveWellWithZeroRate(const Simulator& simulator,
818  const double dt,
819  const GroupStateHelperType& groupStateHelper,
820  WellStateType& well_state)
821  {
822  OPM_TIMEFUNCTION();
823 
824  // Solve a well as stopped with isolation (empty group state for assembly)
825  const auto well_status_orig = this->wellStatus_;
826  this->stopWell();
827 
828  auto inj_controls = Well::InjectionControls(0);
829  auto prod_controls = Well::ProductionControls(0);
830 
831  // Solve with well isolation - the flag "solving_with_zero_rate=true" will be passed down to
832  // assembleWellEqWithoutIterationImpl() to create an empty group state when assembling the
833  // well equations.
834  const bool converged = this->iterateWellEqWithSwitching(
835  simulator, dt, inj_controls, prod_controls,
836  groupStateHelper,
837  well_state,
838  /*fixed_control*/true,
839  /*fixed_status*/true,
840  /*solving_with_zero_rate*/true
841  );
842  this->wellStatus_ = well_status_orig;
843  return converged;
844  }
845 
846  template<typename TypeTag>
847  bool
848  WellInterface<TypeTag>::
849  solveWellForTesting(const Simulator& simulator,
850  const GroupStateHelperType& groupStateHelper,
851  WellStateType& well_state)
852  {
853  OPM_TIMEFUNCTION();
854  auto& deferred_logger = groupStateHelper.deferredLogger();
855 
856  const double dt = simulator.timeStepSize();
857 
858  const auto& summary_state = simulator.vanguard().summaryState();
859  auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
860  auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
861  this->onlyKeepBHPandTHPcontrols(summary_state, well_state, inj_controls, prod_controls);
862 
863  bool converged = false;
864  try {
865  // TODO: the following two functions will be refactored to be one to reduce the code duplication
866  if (!this->param_.local_well_solver_control_switching_){
867  converged = this->iterateWellEqWithControl(
868  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state
869  );
870  } else {
871  if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.well(this->index_of_well_).status == WellStatus::OPEN)) {
872  converged = this->solveWellWithOperabilityCheck(
873  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state
874  );
875  } else {
876  converged = this->iterateWellEqWithSwitching(
877  simulator, dt, inj_controls, prod_controls, groupStateHelper, well_state,
878  /*fixed_control=*/false,
879  /*fixed_status=*/false,
880  /*solving_with_zero_rate=*/false
881  );
882  }
883  }
884 
885  } catch (NumericalProblem& e ) {
886  const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
887  deferred_logger.warning("INNER_ITERATION_FAILED", msg);
888  converged = false;
889  }
890 
891  if (converged) {
892  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " converged");
893  return true;
894  }
895  const int max_iter = this->param_.max_welleq_iter_;
896  deferred_logger.debug("WellTest: Well equation for well " + this->name() + " failed converging in "
897  + std::to_string(max_iter) + " iterations");
898  return false;
899  }
900 
901 
902  template<typename TypeTag>
903  void
904  WellInterface<TypeTag>::
905  solveWellEquation(const Simulator& simulator,
906  const GroupStateHelperType& groupStateHelper,
907  WellStateType& well_state)
908  {
909  OPM_TIMEFUNCTION();
910  auto& deferred_logger = groupStateHelper.deferredLogger();
911  if (!this->isOperableAndSolvable() && !this->wellIsStopped())
912  return;
913 
914  // keep a copy of the original well state
915  const WellStateType well_state0 = well_state;
916  const double dt = simulator.timeStepSize();
917  bool converged = iterateWellEquations(simulator, dt, groupStateHelper, well_state);
918 
919  // Newly opened wells with THP control sometimes struggles to
920  // converge due to bad initial guess. Or due to the simple fact
921  // that the well needs to change to another control.
922  // We therefore try to solve the well with BHP control to get
923  // an better initial guess.
924  // If the well is supposed to operate under THP control
925  // "updateWellControl" will switch it back to THP later.
926  if (!converged) {
927  auto& ws = well_state.well(this->indexOfWell());
928  bool thp_control = false;
929  if (this->well_ecl_.isInjector()) {
930  thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
931  if (thp_control) {
932  ws.injection_cmode = Well::InjectorCMode::BHP;
933  if (this->well_control_log_.empty()) { // only log the first control
934  this->well_control_log_.push_back(WellInjectorCMode2String(Well::InjectorCMode::THP));
935  }
936  }
937  } else {
938  thp_control = ws.production_cmode == Well::ProducerCMode::THP;
939  if (thp_control) {
940  ws.production_cmode = Well::ProducerCMode::BHP;
941  if (this->well_control_log_.empty()) { // only log the first control
942  this->well_control_log_.push_back(WellProducerCMode2String(Well::ProducerCMode::THP));
943  }
944  }
945  }
946  if (thp_control) {
947  const std::string msg = std::string("The newly opened well ") + this->name()
948  + std::string(" with THP control did not converge during inner iterations, we try again with bhp control");
949  deferred_logger.debug(msg);
950  converged = this->iterateWellEquations(simulator, dt, groupStateHelper, well_state);
951  }
952  }
953 
954  if (!converged) {
955  const int max_iter = this->param_.max_welleq_iter_;
956  deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
957  + std::to_string(max_iter) + " iterations");
958  well_state = well_state0;
959  }
960  }
961 
962 
963 
964  template <typename TypeTag>
965  void
966  WellInterface<TypeTag>::
967  assembleWellEq(const Simulator& simulator,
968  const double dt,
969  const GroupStateHelperType& groupStateHelper,
970  WellStateType& well_state)
971  {
972  OPM_TIMEFUNCTION();
973  prepareWellBeforeAssembling(simulator, dt, groupStateHelper, well_state);
974  assembleWellEqWithoutIteration(simulator, groupStateHelper, dt, well_state,
975  /*solving_with_zero_rate=*/false);
976  }
977 
978 
979 
980  template <typename TypeTag>
981  void
982  WellInterface<TypeTag>::
983  assembleWellEqWithoutIteration(const Simulator& simulator,
984  const GroupStateHelperType& groupStateHelper,
985  const double dt,
986  WellStateType& well_state,
987  const bool solving_with_zero_rate)
988  {
989  OPM_TIMEFUNCTION();
990  const auto& summary_state = simulator.vanguard().summaryState();
991  const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
992  const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
993  // TODO: the reason to have inj_controls and prod_controls in the arguments, is that we want to change the control used for the well functions
994  // TODO: maybe we can use std::optional or pointers to simplify here
995  assembleWellEqWithoutIteration(simulator, groupStateHelper, dt, inj_controls, prod_controls, well_state, solving_with_zero_rate);
996  }
997 
998 
999 
1000  template<typename TypeTag>
1001  void
1002  WellInterface<TypeTag>::
1003  updateGroupTargetFallbackFlag(WellStateType& well_state,
1004  DeferredLogger& deferred_logger) const
1005  {
1006  // Get scaled well fractions from derived class
1007  std::vector<Scalar> scaled_well_fractions(FluidSystem::numPhases, 0.0);
1008  this->getScaledWellFractions(scaled_well_fractions, deferred_logger);
1009  // Call the base class method
1010  this->Base::updateGroupTargetFallbackFlag(well_state, scaled_well_fractions, deferred_logger);
1011  }
1012 
1013 
1014 
1015  template<typename TypeTag>
1016  void
1017  WellInterface<TypeTag>::
1018  prepareWellBeforeAssembling(const Simulator& simulator,
1019  const double dt,
1020  const GroupStateHelperType& groupStateHelper,
1021  WellStateType& well_state)
1022  {
1023  OPM_TIMEFUNCTION();
1024  auto& deferred_logger = groupStateHelper.deferredLogger();
1025  const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
1026 
1027  if (this->param_.check_well_operability_iter_)
1028  checkWellOperability(simulator, well_state, groupStateHelper);
1029 
1030  // only use inner well iterations for the first newton iterations.
1031  const auto& iterCtx = simulator.problem().iterationContext();
1032  if (iterCtx.shouldRunInnerWellIterations(this->param_.max_niter_inner_well_iter_)) {
1033  const auto& ws = well_state.well(this->indexOfWell());
1034  const bool nonzero_rate_original =
1035  std::any_of(ws.surface_rates.begin(),
1036  ws.surface_rates.begin() + well_state.numPhases(),
1037  [](Scalar rate) { return rate != Scalar(0.0); });
1038 
1039  this->operability_status_.solvable = true;
1040  if (number_of_well_reopenings_ >= this->param_.max_well_status_switch_) {
1041  // only output the first time
1042  if (number_of_well_reopenings_ == this->param_.max_well_status_switch_) {
1043  const std::string msg = fmt::format("well {} is oscillating between open and stop. \n"
1044  "We don't allow for more than {} re-openings "
1045  "and the well is therefore kept stopped.",
1046  this->name(), number_of_well_reopenings_);
1047  deferred_logger.debug(msg);
1048  changed_to_stopped_this_step_ = old_well_operable;
1049  } else {
1050  changed_to_stopped_this_step_ = false;
1051  }
1052  this->stopWell();
1053  bool converged_zero_rate = this->solveWellWithZeroRate(
1054  simulator, dt, groupStateHelper, well_state
1055  );
1056  if (this->param_.shut_unsolvable_wells_ && !converged_zero_rate ) {
1057  this->operability_status_.solvable = false;
1058  } else {
1059  this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1060  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1061  }
1062  // we increse the number of reopenings to avoid output in the next iteration
1063  number_of_well_reopenings_++;
1064  return;
1065  }
1066  bool converged = this->iterateWellEquations(
1067  simulator, dt, groupStateHelper, well_state
1068  );
1069 
1070  if (converged) {
1071  const bool zero_target = this->wellUnderZeroRateTarget(groupStateHelper);
1072  if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
1073  // Well had non-zero rate, but was stopped during local well-solve. We re-open the well
1074  // for the next global iteration, but if the zero rate persists, it will be stopped.
1075  // This logic is introduced to prevent/ameliorate stopped/revived oscillations
1076  this->operability_status_.resetOperability();
1077  this->openWell();
1078  deferred_logger.debug(" " + this->name() + " is re-opened after being stopped during local solve");
1079  number_of_well_reopenings_++;
1080  }
1081  } else {
1082  // unsolvable wells are treated as not operable and will not be solved for in this iteration.
1083  if (this->param_.shut_unsolvable_wells_) {
1084  this->operability_status_.solvable = false;
1085  }
1086  }
1087  }
1088  if (this->operability_status_.has_negative_potentials) {
1089  auto well_state_copy = well_state;
1090  std::vector<Scalar> potentials;
1091  try {
1092  computeWellPotentials(simulator, well_state_copy, groupStateHelper, potentials);
1093  } catch (const std::exception& e) {
1094  const std::string msg = fmt::format("well {}: computeWellPotentials() failed "
1095  "during attempt to recompute potentials for well: ",
1096  this->name(), e.what());
1097  deferred_logger.info(msg);
1098  this->operability_status_.has_negative_potentials = true;
1099  }
1100  auto& ws = well_state.well(this->indexOfWell());
1101  const int np = well_state.numPhases();
1102  for (int p = 0; p < np; ++p) {
1103  ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
1104  }
1105  }
1106  this->changed_to_open_this_step_ = false;
1107  changed_to_stopped_this_step_ = false;
1108 
1109  const bool well_operable = this->operability_status_.isOperableAndSolvable();
1110  if (!well_operable) {
1111  this->stopWell();
1112  try {
1113  this->solveWellWithZeroRate(
1114  simulator, dt, groupStateHelper, well_state
1115  );
1116  } catch (const std::exception& e) {
1117  const std::string msg = fmt::format("well {}: solveWellWithZeroRate() failed "
1118  "during attempt to solve with zero rate for well: ",
1119  this->name(), e.what());
1120  deferred_logger.info(msg);
1121  // we set the rate to zero to make sure the well dont contribute to the group rate
1122  auto& ws = well_state.well(this->indexOfWell());
1123  const int np = well_state.numPhases();
1124  for (int p = 0; p < np; ++p) {
1125  ws.surface_rates[p] = Scalar{0.0};
1126  }
1127  }
1128  if (old_well_operable) {
1129  const std::string ctx = iterCtx.inLocalSolve() ? " (NLDD domain solve)" : "";
1130  deferred_logger.debug(" well " + this->name() + " gets STOPPED during iteration" + ctx);
1131  changed_to_stopped_this_step_ = true;
1132  }
1133  } else if (well_state.isOpen(this->name())) {
1134  this->openWell();
1135  if (!old_well_operable) {
1136  const std::string ctx = iterCtx.inLocalSolve() ? " (NLDD domain solve)" : "";
1137  deferred_logger.debug(" well " + this->name() + " gets REVIVED during iteration" + ctx);
1138  this->changed_to_open_this_step_ = true;
1139  }
1140  }
1141  }
1142 
1143  template<typename TypeTag>
1144  void
1145  WellInterface<TypeTag>::addCellRates(std::map<int, RateVector>& cellRates_) const
1146  {
1147  if(!this->operability_status_.solvable)
1148  return;
1149 
1150  for (int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
1151  const auto cellIdx = this->cells()[perfIdx];
1152  const auto it = cellRates_.find(cellIdx);
1153  RateVector rates = (it == cellRates_.end()) ? 0.0 : it->second;
1154  for (auto i=0*RateVector::dimension; i < RateVector::dimension; ++i)
1155  {
1156  rates[i] += connectionRates_[perfIdx][i];
1157  }
1158  cellRates_.insert_or_assign(cellIdx, rates);
1159  }
1160  }
1161 
1162  template<typename TypeTag>
1163  typename WellInterface<TypeTag>::Scalar
1164  WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const
1165  {
1166  for (int perfIdx = 0; perfIdx < this->number_of_local_perforations_; ++perfIdx) {
1167  if (this->cells()[perfIdx] == cellIdx) {
1168  const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
1169  return connectionRates_[perfIdx][activeCompIdx].value();
1170  }
1171  }
1172  // this is not thread safe
1173  OPM_THROW(std::invalid_argument, "The well with name " + this->name()
1174  + " does not perforate cell " + std::to_string(cellIdx));
1175  return 0.0;
1176  }
1177 
1178 
1179 
1180 
1181  template<typename TypeTag>
1182  void
1183  WellInterface<TypeTag>::
1184  checkWellOperability(const Simulator& simulator,
1185  const WellStateType& well_state,
1186  const GroupStateHelperType& groupStateHelper)
1187  {
1188  auto& deferred_logger = groupStateHelper.deferredLogger();
1189  OPM_TIMEFUNCTION();
1190  if (!this->param_.check_well_operability_) {
1191  return;
1192  }
1193 
1194  if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
1195  return;
1196  }
1197 
1198  updateWellOperability(simulator, well_state, groupStateHelper);
1199  if (!this->operability_status_.isOperableAndSolvable()) {
1200  this->operability_status_.use_vfpexplicit = true;
1201  deferred_logger.debug("EXPLICIT_LOOKUP_VFP",
1202  "well not operable, trying with explicit vfp lookup: " + this->name());
1203  updateWellOperability(simulator, well_state, groupStateHelper);
1204  }
1205  }
1206 
1207 
1208 
1209  template<typename TypeTag>
1210  void
1211  WellInterface<TypeTag>::
1212  gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator,
1213  WellStateType& well_state,
1214  const GroupState<Scalar>& group_state,
1215  GLiftEclWells& ecl_well_map,
1216  DeferredLogger& deferred_logger)
1217  {
1218  OPM_TIMEFUNCTION();
1219  const auto& summary_state = simulator.vanguard().summaryState();
1220  const auto& well_name = this->name();
1221  if (!this->wellHasTHPConstraints(summary_state)) {
1222  const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
1223  deferred_logger.info(msg);
1224  return;
1225  }
1226  const auto& schedule = simulator.vanguard().schedule();
1227  const auto report_step_idx = simulator.episodeIndex();
1228  const auto& glo = schedule.glo(report_step_idx);
1229  if (!glo.has_well(well_name)) {
1230  const std::string msg = fmt::format(
1231  "GLIFT WTEST: Well {} : Gas lift not activated: "
1232  "WLIFTOPT is probably missing. Skipping.", well_name);
1233  deferred_logger.info(msg);
1234  return;
1235  }
1236  const auto& gl_well = glo.well(well_name);
1237 
1238  // Use gas lift optimization to get ALQ for well test
1239  std::unique_ptr<GasLiftSingleWell> glift =
1240  initializeGliftWellTest_<GasLiftSingleWell>(simulator,
1241  well_state,
1242  group_state,
1243  ecl_well_map,
1244  deferred_logger);
1245  auto [wtest_alq, success] = glift->wellTestALQ();
1246  std::string msg;
1247  const auto& unit_system = schedule.getUnits();
1248  if (success) {
1249  well_state.well(well_name).alq_state.set(wtest_alq);
1250  msg = fmt::format(
1251  "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
1252  well_name, unit_system.from_si(UnitSystem::measure::gas_surface_rate, wtest_alq));
1253  }
1254  else {
1255  if (!gl_well.use_glo()) {
1256  msg = fmt::format(
1257  "GLIFT WTEST: Well {} : Gas lift optimization deactivated. Setting ALQ to WLIFTOPT item 3 = {}",
1258  well_name,
1259  unit_system.from_si(UnitSystem::measure::gas_surface_rate, well_state.well(well_name).alq_state.get()));
1260 
1261  }
1262  else {
1263  msg = fmt::format(
1264  "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
1265  well_name);
1266  }
1267  }
1268  deferred_logger.info(msg);
1269  }
1270 
1271  template<typename TypeTag>
1272  void
1273  WellInterface<TypeTag>::
1274  updateWellOperability(const Simulator& simulator,
1275  const WellStateType& well_state,
1276  const GroupStateHelperType& groupStateHelper)
1277  {
1278  auto& deferred_logger = groupStateHelper.deferredLogger();
1279  OPM_TIMEFUNCTION();
1280  if (this->param_.local_well_solver_control_switching_) {
1281  const bool success = updateWellOperabilityFromWellEq(simulator, groupStateHelper);
1282  if (!success) {
1283  this->operability_status_.solvable = false;
1284  deferred_logger.debug("Operability check using well equations did not converge for well "
1285  + this->name() + ". Mark the well as unsolvable." );
1286  }
1287  return;
1288  }
1289  this->operability_status_.resetOperability();
1290 
1291  bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
1292  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
1293  bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
1294  well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
1295 
1296  // Operability checking is not free
1297  // Only check wells under BHP and THP control
1298  bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
1299  if (check_thp || bhp_controlled) {
1300  updateIPR(simulator, deferred_logger);
1301  checkOperabilityUnderBHPLimit(well_state, simulator, deferred_logger);
1302  }
1303  // we do some extra checking for wells under THP control.
1304  if (check_thp) {
1305  checkOperabilityUnderTHPLimit(simulator, well_state, groupStateHelper);
1306  }
1307  }
1308 
1309  template<typename TypeTag>
1310  bool
1311  WellInterface<TypeTag>::
1312  updateWellOperabilityFromWellEq(const Simulator& simulator,
1313  const GroupStateHelperType& groupStateHelper)
1314  {
1315  OPM_TIMEFUNCTION();
1316  // only makes sense if we're using this parameter is true
1317  assert(this->param_.local_well_solver_control_switching_);
1318  this->operability_status_.resetOperability();
1319  GroupStateHelperType groupStateHelper_copy = groupStateHelper;
1320  WellStateType well_state_copy = groupStateHelper_copy.wellState();
1321  const double dt = simulator.timeStepSize();
1322  // Ensure that groupStateHelper uses well_state_copy as WellState for iterateWellEquations()
1323  // and the guard ensures that the original well state is restored at scope exit, i.e. at
1324  // the end of this function.
1325  auto guard = groupStateHelper_copy.pushWellState(well_state_copy);
1326  // equations should be converged at this stage, so only one it is needed
1327  bool converged = iterateWellEquations(simulator, dt, groupStateHelper_copy, well_state_copy);
1328  return converged;
1329  }
1330 
1331  template<typename TypeTag>
1332  void
1333  WellInterface<TypeTag>::
1334  scaleSegmentRatesAndPressure([[maybe_unused]] WellStateType& well_state) const
1335  {
1336  // only relevant for MSW
1337  }
1338 
1339  template<typename TypeTag>
1340  void
1341  WellInterface<TypeTag>::
1342  updateWellStateWithTarget(const Simulator& simulator,
1343  const GroupStateHelperType& groupStateHelper,
1344  WellStateType& well_state) const
1345  {
1346  OPM_TIMEFUNCTION();
1347  auto& deferred_logger = groupStateHelper.deferredLogger();
1348  // only bhp and wellRates are used to initilize the primaryvariables for standard wells
1349  const auto& well = this->well_ecl_;
1350  const int well_index = this->index_of_well_;
1351  auto& ws = well_state.well(well_index);
1352  const int np = well_state.numPhases();
1353  const auto& summaryState = simulator.vanguard().summaryState();
1354  const auto& schedule = simulator.vanguard().schedule();
1355 
1356  // Discard old primary variables, the new well state
1357  // may not be anywhere near the old one.
1358  ws.primaryvar.resize(0);
1359 
1360  if (this->wellIsStopped()) {
1361  for (int p = 0; p<np; ++p) {
1362  ws.surface_rates[p] = 0;
1363  }
1364  ws.thp = 0;
1365  return;
1366  }
1367 
1368  if (this->isInjector() )
1369  {
1370  const auto& controls = well.injectionControls(summaryState);
1371 
1372  InjectorType injectorType = controls.injector_type;
1373  int phasePos;
1374  switch (injectorType) {
1375  case InjectorType::WATER:
1376  {
1377  phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1378  break;
1379  }
1380  case InjectorType::OIL:
1381  {
1382  phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1383  break;
1384  }
1385  case InjectorType::GAS:
1386  {
1387  phasePos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1388  break;
1389  }
1390  default:
1391  OPM_DEFLOG_THROW(std::runtime_error, "Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
1392  }
1393 
1394  const auto current = ws.injection_cmode;
1395 
1396  switch (current) {
1397  case Well::InjectorCMode::RATE:
1398  {
1399  ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
1400  if(this->rsRvInj() > 0) {
1401  if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1402  const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1403  ws.surface_rates[gas_pos] = controls.surface_rate * this->rsRvInj();
1404  } else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1405  const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1406  ws.surface_rates[oil_pos] = controls.surface_rate * this->rsRvInj();
1407  } else {
1408  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 );
1409  }
1410  }
1411  break;
1412  }
1413 
1414  case Well::InjectorCMode::RESV:
1415  {
1416  std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1417  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
1418  const Scalar coeff = convert_coeff[phasePos];
1419  ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
1420  break;
1421  }
1422 
1423  case Well::InjectorCMode::THP:
1424  {
1425  auto rates = ws.surface_rates;
1426  Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
1427  rates,
1428  well,
1429  summaryState,
1430  this->getRefDensity(),
1431  deferred_logger);
1432  ws.bhp = bhp;
1433  ws.thp = this->getTHPConstraint(summaryState);
1434 
1435  // if the total rates are negative or zero
1436  // we try to provide a better intial well rate
1437  // using the well potentials
1438  Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
1439  if (total_rate <= 0.0)
1440  ws.surface_rates = ws.well_potentials;
1441 
1442  break;
1443  }
1444  case Well::InjectorCMode::BHP:
1445  {
1446  ws.bhp = controls.bhp_limit;
1447  Scalar total_rate = 0.0;
1448  for (int p = 0; p<np; ++p) {
1449  total_rate += ws.surface_rates[p];
1450  }
1451  // if the total rates are negative or zero
1452  // we try to provide a better intial well rate
1453  // using the well potentials
1454  if (total_rate <= 0.0)
1455  ws.surface_rates = ws.well_potentials;
1456 
1457  break;
1458  }
1459  case Well::InjectorCMode::GRUP:
1460  {
1461  assert(well.isAvailableForGroupControl());
1462  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1463  const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1464  well_state[well.name()].efficiency_scaling_factor;
1465  std::optional<Scalar> target =
1466  this->getGroupInjectionTargetRate(group,
1467  groupStateHelper,
1468  injectorType,
1469  efficiencyFactor);
1470  if (target)
1471  ws.surface_rates[phasePos] = *target;
1472  break;
1473  }
1474  case Well::InjectorCMode::CMODE_UNDEFINED:
1475  {
1476  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger );
1477  }
1478 
1479  }
1480  // for wells with zero injection rate, if we assign exactly zero rate,
1481  // we will have to assume some trivial composition in the wellbore.
1482  // here, we use some small value (about 0.01 m^3/day ~= 1.e-7) to initialize
1483  // the zero rate target, then we can use to retain the composition information
1484  // within the wellbore from the previous result, and hopefully it is a good
1485  // initial guess for the zero rate target.
1486  ws.surface_rates[phasePos] = std::max(Scalar{1.e-7}, ws.surface_rates[phasePos]);
1487 
1488  if (ws.bhp == 0.) {
1489  ws.bhp = controls.bhp_limit;
1490  }
1491  }
1492  //Producer
1493  else
1494  {
1495  const auto current = ws.production_cmode;
1496  const auto& controls = well.productionControls(summaryState);
1497  switch (current) {
1498  case Well::ProducerCMode::ORAT:
1499  {
1500  const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1501  Scalar current_rate = -ws.surface_rates[oil_pos];
1502  // for trivial rates or opposite direction we don't just scale the rates
1503  // but use either the potentials or the mobility ratio to initial the well rates
1504  if (current_rate > 0.0) {
1505  for (int p = 0; p<np; ++p) {
1506  ws.surface_rates[p] *= controls.oil_rate/current_rate;
1507  }
1508  } else {
1509  const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1510  double control_fraction = fractions[oil_pos];
1511  if (control_fraction != 0.0) {
1512  for (int p = 0; p<np; ++p) {
1513  ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
1514  }
1515  }
1516  }
1517  break;
1518  }
1519  case Well::ProducerCMode::WRAT:
1520  {
1521  const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1522  Scalar current_rate = -ws.surface_rates[water_pos];
1523  // for trivial rates or opposite direction we don't just scale the rates
1524  // but use either the potentials or the mobility ratio to initial the well rates
1525  if (current_rate > 0.0) {
1526  for (int p = 0; p<np; ++p) {
1527  ws.surface_rates[p] *= controls.water_rate/current_rate;
1528  }
1529  } else {
1530  const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1531  const Scalar control_fraction = fractions[water_pos];
1532  if (control_fraction != 0.0) {
1533  for (int p = 0; p<np; ++p) {
1534  ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
1535  }
1536  }
1537  }
1538  break;
1539  }
1540  case Well::ProducerCMode::GRAT:
1541  {
1542  const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1543  Scalar current_rate = -ws.surface_rates[gas_pos];
1544  // or trivial rates or opposite direction we don't just scale the rates
1545  // but use either the potentials or the mobility ratio to initial the well rates
1546  if (current_rate > 0.0) {
1547  for (int p = 0; p<np; ++p) {
1548  ws.surface_rates[p] *= controls.gas_rate/current_rate;
1549  }
1550  } else {
1551  const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
1552  const Scalar control_fraction = fractions[gas_pos];
1553  if (control_fraction != 0.0) {
1554  for (int p = 0; p<np; ++p) {
1555  ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
1556  }
1557  }
1558  }
1559 
1560  break;
1561 
1562  }
1563  case Well::ProducerCMode::LRAT:
1564  {
1565  const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1566  const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1567  Scalar current_rate = - ws.surface_rates[water_pos]
1568  - ws.surface_rates[oil_pos];
1569  // or trivial rates or opposite direction we don't just scale the rates
1570  // but use either the potentials or the mobility ratio to initial the well rates
1571  if (current_rate > 0.0) {
1572  for (int p = 0; p<np; ++p) {
1573  ws.surface_rates[p] *= controls.liquid_rate/current_rate;
1574  }
1575  } else {
1576  const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1577  const Scalar control_fraction = fractions[water_pos] + fractions[oil_pos];
1578  if (control_fraction != 0.0) {
1579  for (int p = 0; p<np; ++p) {
1580  ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
1581  }
1582  }
1583  }
1584  break;
1585  }
1586  case Well::ProducerCMode::CRAT:
1587  {
1588  OPM_DEFLOG_THROW(std::runtime_error,
1589  fmt::format("CRAT control not supported, well {}", this->name()),
1590  deferred_logger);
1591  }
1592  case Well::ProducerCMode::RESV:
1593  {
1594  std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
1595  this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
1596  Scalar total_res_rate = 0.0;
1597  for (int p = 0; p<np; ++p) {
1598  total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
1599  }
1600  if (controls.prediction_mode) {
1601  // or trivial rates or opposite direction we don't just scale the rates
1602  // but use either the potentials or the mobility ratio to initial the well rates
1603  if (total_res_rate > 0.0) {
1604  for (int p = 0; p<np; ++p) {
1605  ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
1606  }
1607  } else {
1608  const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1609  for (int p = 0; p<np; ++p) {
1610  ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
1611  }
1612  }
1613  } else {
1614  std::vector<Scalar> hrates(this->number_of_phases_,0.);
1615  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1616  const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
1617  hrates[phase_pos] = controls.water_rate;
1618  }
1619  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1620  const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
1621  hrates[phase_pos] = controls.oil_rate;
1622  }
1623  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1624  const int phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1625  hrates[phase_pos] = controls.gas_rate;
1626  }
1627  std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
1628  this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
1629  Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1630  // or trivial rates or opposite direction we don't just scale the rates
1631  // but use either the potentials or the mobility ratio to initial the well rates
1632  if (total_res_rate > 0.0) {
1633  for (int p = 0; p<np; ++p) {
1634  ws.surface_rates[p] *= target/total_res_rate;
1635  }
1636  } else {
1637  const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
1638  for (int p = 0; p<np; ++p) {
1639  ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1640  }
1641  }
1642  }
1643  break;
1644  }
1645  case Well::ProducerCMode::BHP:
1646  {
1647  ws.bhp = controls.bhp_limit;
1648  Scalar total_rate = 0.0;
1649  for (int p = 0; p<np; ++p) {
1650  total_rate -= ws.surface_rates[p];
1651  }
1652  // if the total rates are negative or zero
1653  // we try to provide a better intial well rate
1654  // using the well potentials
1655  if (total_rate <= 0.0){
1656  for (int p = 0; p<np; ++p) {
1657  ws.surface_rates[p] = -ws.well_potentials[p];
1658  }
1659  }
1660  break;
1661  }
1662  case Well::ProducerCMode::THP:
1663  {
1664  const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, groupStateHelper);
1665 
1666  if (!update_success) {
1667  // the following is the original way of initializing well state with THP constraint
1668  // keeping it for robust reason in case that it fails to get a bhp value with THP constraint
1669  // more sophisticated design might be needed in the future
1670  auto rates = ws.surface_rates;
1671  this->adaptRatesForVFP(rates);
1672  const Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
1673  well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1674  ws.bhp = bhp;
1675  ws.thp = this->getTHPConstraint(summaryState);
1676  // if the total rates are negative or zero
1677  // we try to provide a better initial well rate
1678  // using the well potentials
1679  const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1680  if (total_rate <= 0.0) {
1681  for (int p = 0; p < this->number_of_phases_; ++p) {
1682  ws.surface_rates[p] = -ws.well_potentials[p];
1683  }
1684  }
1685  }
1686  break;
1687  }
1688  case Well::ProducerCMode::GRUP:
1689  {
1690  assert(well.isAvailableForGroupControl());
1691  this->updateGroupTargetFallbackFlag(well_state, deferred_logger);
1692  const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1693  const Scalar efficiencyFactor = well.getEfficiencyFactor() *
1694  well_state[well.name()].efficiency_scaling_factor;
1695  Scalar scale = this->getGroupProductionTargetRate(group,
1696  groupStateHelper,
1697  efficiencyFactor);
1698 
1699  // we don't want to scale with zero and get zero rates.
1700  if (scale > 0) {
1701  for (int p = 0; p<np; ++p) {
1702  ws.surface_rates[p] *= scale;
1703  }
1704  ws.trivial_group_target = false;
1705  } else {
1706  // If group target is trivial we dont want to flip to other controls. To avoid oscillation we store
1707  // this information in the well state and explicitly check for this condition when evaluating well controls.
1708  ws.trivial_group_target = true;
1709  }
1710  break;
1711  }
1712  case Well::ProducerCMode::CMODE_UNDEFINED:
1713  case Well::ProducerCMode::NONE:
1714  {
1715  OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name() , deferred_logger);
1716  break;
1717  }
1718  } // end of switch
1719 
1720  if (ws.bhp == 0.) {
1721  ws.bhp = controls.bhp_limit;
1722  }
1723  }
1724  }
1725 
1726  template<typename TypeTag>
1727  bool
1728  WellInterface<TypeTag>::
1729  wellUnderZeroRateTarget(const GroupStateHelperType& groupStateHelper) const
1730  {
1731  OPM_TIMEFUNCTION();
1732  const auto& well_state = groupStateHelper.wellState();
1733  // Check if well is under zero rate control, either directly or from group
1734  const bool isGroupControlled = this->wellUnderGroupControl(well_state.well(this->index_of_well_));
1735  if (!isGroupControlled) {
1736  // well is not under group control, check "individual" version
1737  const auto& summaryState = groupStateHelper.summaryState();
1738  return this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
1739  } else {
1740  return this->wellUnderZeroGroupRateTarget(groupStateHelper, isGroupControlled);
1741  }
1742  }
1743 
1744  template <typename TypeTag>
1745  bool
1746  WellInterface<TypeTag>::wellUnderZeroGroupRateTarget(const GroupStateHelperType& groupStateHelper,
1747  const std::optional<bool> group_control) const
1748  {
1749  const auto& well_state = groupStateHelper.wellState();
1750  // Check if well is under zero rate target from group
1751  const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.well(this->index_of_well_)));
1752  if (isGroupControlled) {
1753  return this->zeroGroupRateTarget(groupStateHelper);
1754  }
1755  return false;
1756  }
1757 
1758  template<typename TypeTag>
1759  bool
1760  WellInterface<TypeTag>::
1761  stoppedOrZeroRateTarget(const GroupStateHelperType& groupStateHelper) const
1762  {
1763  // Check if well is stopped or under zero rate control, either
1764  // directly or from group.
1765  return this->wellIsStopped()
1766  || this->wellUnderZeroRateTarget(groupStateHelper);
1767  }
1768 
1769  template<typename TypeTag>
1770  std::vector<typename WellInterface<TypeTag>::Scalar>
1771  WellInterface<TypeTag>::
1772  initialWellRateFractions(const Simulator& simulator,
1773  const WellStateType& well_state) const
1774  {
1775  OPM_TIMEFUNCTION();
1776  const int np = this->number_of_phases_;
1777  std::vector<Scalar> scaling_factor(np);
1778  const auto& ws = well_state.well(this->index_of_well_);
1779 
1780  Scalar total_potentials = 0.0;
1781  for (int p = 0; p<np; ++p) {
1782  total_potentials += ws.well_potentials[p];
1783  }
1784  if (total_potentials > 0) {
1785  for (int p = 0; p<np; ++p) {
1786  scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1787  }
1788  return scaling_factor;
1789  }
1790  // if we don't have any potentials we weight it using the mobilites
1791  // We only need approximation so we don't bother with the vapporized oil and dissolved gas
1792  Scalar total_tw = 0;
1793  const int nperf = this->number_of_local_perforations_;
1794  for (int perf = 0; perf < nperf; ++perf) {
1795  total_tw += this->well_index_[perf];
1796  }
1797  total_tw = this->parallelWellInfo().communication().sum(total_tw);
1798 
1799  for (int perf = 0; perf < nperf; ++perf) {
1800  const int cell_idx = this->well_cells_[perf];
1801  const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1802  const auto& fs = intQuants.fluidState();
1803  const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
1804  Scalar total_mobility = 0.0;
1805  for (int p = 0; p < np; ++p) {
1806  const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1807  total_mobility += fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value();
1808  }
1809  for (int p = 0; p < np; ++p) {
1810  const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
1811  scaling_factor[p] += well_tw_fraction * fs.invB(canonical_phase_idx).value() * intQuants.mobility(canonical_phase_idx).value() / total_mobility;
1812  }
1813  }
1814  return scaling_factor;
1815  }
1816 
1817 
1818 
1819  template <typename TypeTag>
1820  void
1821  WellInterface<TypeTag>::
1822  initializeProducerWellState(const Simulator& simulator,
1823  WellStateType& well_state,
1824  DeferredLogger& deferred_logger) const
1825  {
1826  assert(this->isProducer());
1827  OPM_TIMEFUNCTION();
1828  // Check if the rates of this well only are single-phase, do nothing
1829  // if more than one nonzero rate.
1830  auto& ws = well_state.well(this->index_of_well_);
1831  int nonzero_rate_index = -1;
1832  const Scalar floating_point_error_epsilon = 1e-14;
1833  for (int p = 0; p < this->number_of_phases_; ++p) {
1834  if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1835  if (nonzero_rate_index == -1) {
1836  nonzero_rate_index = p;
1837  } else {
1838  // More than one nonzero rate.
1839  return;
1840  }
1841  }
1842  }
1843 
1844  // Calculate rates at bhp limit, or 1 bar if no limit.
1845  std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
1846  bool rates_evaluated_at_1bar = false;
1847  {
1848  const auto& summary_state = simulator.vanguard().summaryState();
1849  const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
1850  const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
1851  this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
1852  // Remember of we evaluated the rates at (approx.) 1 bar or not.
1853  rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
1854  // Check that no rates are positive.
1855  if (std::ranges::any_of(well_q_s, [](Scalar q) { return q > 0.0; })) {
1856  // Did we evaluate at 1 bar? If not, then we can try again at 1 bar.
1857  if (!rates_evaluated_at_1bar) {
1858  this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
1859  rates_evaluated_at_1bar = true;
1860  }
1861  // At this point we can only set the wrong-direction (if any) values to zero.
1862  for (auto& q : well_q_s) {
1863  q = std::min(q, Scalar{0.0});
1864  }
1865  }
1866  }
1867 
1868  if (nonzero_rate_index == -1) {
1869  // No nonzero rates on input.
1870  // Use the computed rate directly, or scaled by a factor
1871  // 0.5 (to avoid too high values) if it was evaluated at 1 bar.
1872  const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
1873  for (int p = 0; p < this->number_of_phases_; ++p) {
1874  ws.surface_rates[p] = factor * well_q_s[p];
1875  }
1876  return;
1877  }
1878 
1879  // If we are here, we had a single nonzero rate for the well,
1880  // typically from a rate constraint. We must make sure it is
1881  // respected, so if it was lower than the calculated rate for
1882  // the same phase we scale all rates to match.
1883  const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1884  const Scalar computed_rate = well_q_s[nonzero_rate_index];
1885  if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {
1886  // Note that both rates below are negative. The factor should be < 1.0.
1887  const Scalar factor = initial_nonzero_rate / computed_rate;
1888  assert(factor < 1.0);
1889  for (int p = 0; p < this->number_of_phases_; ++p) {
1890  // We skip the nonzero_rate_index, as that should remain as it was.
1891  if (p != nonzero_rate_index) {
1892  ws.surface_rates[p] = factor * well_q_s[p];
1893  }
1894  }
1895  return;
1896  }
1897 
1898  // If we are here, we had a single nonzero rate, but it was
1899  // higher than the one calculated from the bhp limit, so we
1900  // use the calculated rates.
1901  for (int p = 0; p < this->number_of_phases_; ++p) {
1902  ws.surface_rates[p] = well_q_s[p];
1903  }
1904  }
1905 
1906  template <typename TypeTag>
1907  template<class Value>
1908  void
1909  WellInterface<TypeTag>::
1910  getTw(std::vector<Value>& Tw,
1911  const int perf,
1912  const IntensiveQuantities& intQuants,
1913  const Value& trans_mult,
1914  const SingleWellStateType& ws) const
1915  {
1916  OPM_TIMEFUNCTION_LOCAL(Subsystem::Wells);
1917  // Add a Forchheimer term to the gas phase CTF if the run uses
1918  // either of the WDFAC or the WDFACCOR keywords.
1919  if (static_cast<std::size_t>(perf) >= this->well_cells_.size()) {
1920  OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly wellIndex was called with a global instead of a local perforation index!");
1921  }
1922 
1923  if constexpr (! Indices::gasEnabled) {
1924  return;
1925  }
1926 
1927  const auto& wdfac = this->well_ecl_.getWDFAC();
1928 
1929  if (! wdfac.useDFactor() || (this->well_index_[perf] == 0.0)) {
1930  return;
1931  }
1932 
1933  const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
1934  if (d < 1.0e-15) {
1935  return;
1936  }
1937 
1938  // Solve quadratic equations for connection rates satisfying the ipr and the flow-dependent skin.
1939  // If more than one solution, pick the one corresponding to lowest absolute rate (smallest skin).
1940  const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
1941  const Scalar Kh = connection.Kh();
1942  const Scalar scaling = std::numbers::pi * Kh * connection.wpimult();
1943  const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
1944 
1945  const Scalar connection_pressure = ws.perf_data.pressure[perf];
1946  const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
1947  const Scalar drawdown = cell_pressure - connection_pressure;
1948  const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
1949  const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
1950  const Scalar a = d;
1951  const Scalar b = 2 * scaling / getValue(Tw[gas_comp_idx]);
1952  const Scalar c = -2 * scaling * mob_g * drawdown;
1953 
1954  Scalar consistent_Q = -1.0e20;
1955  // Find and check negative solutions (a --> -a)
1956  const Scalar r2n = b*b + 4*a*c;
1957  if (r2n >= 0) {
1958  const Scalar rn = std::sqrt(r2n);
1959  const Scalar xn1 = (b-rn)*0.5/a;
1960  if (xn1 <= 0) {
1961  consistent_Q = xn1;
1962  }
1963  const Scalar xn2 = (b+rn)*0.5/a;
1964  if (xn2 <= 0 && xn2 > consistent_Q) {
1965  consistent_Q = xn2;
1966  }
1967  }
1968  // Find and check positive solutions
1969  consistent_Q *= -1;
1970  const Scalar r2p = b*b - 4*a*c;
1971  if (r2p >= 0) {
1972  const Scalar rp = std::sqrt(r2p);
1973  const Scalar xp1 = (rp-b)*0.5/a;
1974  if (xp1 > 0 && xp1 < consistent_Q) {
1975  consistent_Q = xp1;
1976  }
1977  const Scalar xp2 = -(rp+b)*0.5/a;
1978  if (xp2 > 0 && xp2 < consistent_Q) {
1979  consistent_Q = xp2;
1980  }
1981  }
1982  Tw[gas_comp_idx] = 1.0 / (1.0 / (trans_mult * this->well_index_[perf]) + (consistent_Q/2 * d / scaling));
1983  }
1984 
1985  template <typename TypeTag>
1986  void
1987  WellInterface<TypeTag>::
1988  updateConnectionDFactor(const Simulator& simulator,
1989  SingleWellStateType& ws) const
1990  {
1991  if (! this->well_ecl_.getWDFAC().useDFactor()) {
1992  return;
1993  }
1994 
1995  auto& d_factor = ws.perf_data.connection_d_factor;
1996 
1997  for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1998  const int cell_idx = this->well_cells_[perf];
1999  const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2000 
2001  d_factor[perf] = this->computeConnectionDFactor(perf, intQuants, ws);
2002  }
2003  }
2004 
2005  template <typename TypeTag>
2006  typename WellInterface<TypeTag>::Scalar
2007  WellInterface<TypeTag>::
2008  computeConnectionDFactor(const int perf,
2009  const IntensiveQuantities& intQuants,
2010  const SingleWellStateType& ws) const
2011  {
2012  auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
2013  return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
2014  };
2015 
2016  // Viscosity is evaluated at connection pressure.
2017  auto gas_visc = [connection_pressure = ws.perf_data.pressure[perf],
2018  temperature = ws.temperature,
2019  regIdx = this->pvtRegionIdx(), &intQuants]()
2020  {
2021  const auto rv = getValue(intQuants.fluidState().Rv());
2022 
2023  const auto& gasPvt = FluidSystem::gasPvt();
2024 
2025  // Note that rv here is from grid block with typically
2026  // p_block > connection_pressure
2027  // so we may very well have rv > rv_sat
2028  const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
2029  (regIdx, temperature, connection_pressure);
2030 
2031  if (! (rv < rv_sat)) {
2032  return gasPvt.saturatedViscosity(regIdx, temperature,
2033  connection_pressure);
2034  }
2035 
2036  return gasPvt.viscosity(regIdx, temperature, connection_pressure,
2037  rv, getValue(intQuants.fluidState().Rvw()));
2038  };
2039 
2040  const auto& connection = this->well_ecl_.getConnections()
2041  [ws.perf_data.ecl_index[perf]];
2042 
2043  return this->well_ecl_.getWDFAC().getDFactor(rhoGS, gas_visc, connection);
2044  }
2045 
2046 
2047  template <typename TypeTag>
2048  void
2049  WellInterface<TypeTag>::
2050  updateConnectionTransmissibilityFactor(const Simulator& simulator,
2051  SingleWellStateType& ws) const
2052  {
2053  auto connCF = [&connIx = std::as_const(ws.perf_data.ecl_index),
2054  &conns = this->well_ecl_.getConnections()]
2055  (const int perf)
2056  {
2057  return conns[connIx[perf]].CF();
2058  };
2059 
2060  auto obtain = [](const Eval& value)
2061  {
2062  return getValue(value);
2063  };
2064 
2065  auto& tmult = ws.perf_data.connection_compaction_tmult;
2066  auto& ctf = ws.perf_data.connection_transmissibility_factor;
2067 
2068  for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2069  const int cell_idx = this->well_cells_[perf];
2070  Scalar trans_mult(0.0);
2071  getTransMult(trans_mult, simulator, cell_idx, obtain);
2072  tmult[perf] = trans_mult;
2073 
2074  ctf[perf] = connCF(perf) * tmult[perf];
2075  }
2076  }
2077 
2078 
2079  template<typename TypeTag>
2080  typename WellInterface<TypeTag>::Eval
2081  WellInterface<TypeTag>::getPerfCellPressure(const typename WellInterface<TypeTag>::FluidState& fs) const
2082  {
2083  if constexpr (Indices::oilEnabled) {
2084  return fs.pressure(FluidSystem::oilPhaseIdx);
2085  } else if constexpr (Indices::gasEnabled) {
2086  return fs.pressure(FluidSystem::gasPhaseIdx);
2087  } else {
2088  return fs.pressure(FluidSystem::waterPhaseIdx);
2089  }
2090  }
2091 
2092  template <typename TypeTag>
2093  template<class Value, class Callback>
2094  void
2095  WellInterface<TypeTag>::
2096  getTransMult(Value& trans_mult,
2097  const Simulator& simulator,
2098  const int cell_idx,
2099  Callback& extendEval) const
2100  {
2101  const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2102  trans_mult = simulator.problem().template wellTransMultiplier<Value>(intQuants, cell_idx, extendEval);
2103  }
2104 
2105  template <typename TypeTag>
2106  template<class Value, class Callback>
2107  void
2108  WellInterface<TypeTag>::
2109  getMobility(const Simulator& simulator,
2110  const int local_perf_index,
2111  std::vector<Value>& mob,
2112  Callback& extendEval,
2113  [[maybe_unused]] DeferredLogger& deferred_logger) const
2114  {
2115  auto relpermArray = []()
2116  {
2117  if constexpr (std::is_same_v<Value, Scalar>) {
2118  return std::array<Scalar,3>{};
2119  } else {
2120  return std::array<Eval,3>{};
2121  }
2122  };
2123  if (static_cast<std::size_t>(local_perf_index) >= this->well_cells_.size()) {
2124  OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly getMobility was called with a global instead of a local perforation index!");
2125  }
2126  const int cell_idx = this->well_cells_[local_perf_index];
2127  assert (int(mob.size()) == this->num_conservation_quantities_);
2128  const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
2129  const auto& materialLawManager = simulator.problem().materialLawManager();
2130 
2131  // either use mobility of the perforation cell or calculate its own
2132  // based on passing the saturation table index
2133  const int satid = this->saturation_table_number_[local_perf_index] - 1;
2134  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
2135  if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell
2136  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2137  if (!FluidSystem::phaseIsActive(phaseIdx)) {
2138  continue;
2139  }
2140 
2141  const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
2142  mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
2143  }
2144  if constexpr (has_solvent) {
2145  mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility());
2146  }
2147  } else {
2148  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
2149  auto relativePerms = relpermArray();
2150  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
2151 
2152  // reset the satnumvalue back to original
2153  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
2154 
2155  // compute the mobility
2156  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2157  if (!FluidSystem::phaseIsActive(phaseIdx)) {
2158  continue;
2159  }
2160 
2161  const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
2162  mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
2163  }
2164 
2165  if constexpr (has_solvent) {
2166  const auto Fsolgas = intQuants.solventSaturation() / (intQuants.solventSaturation() + intQuants.fluidState().saturation(FluidSystem::gasPhaseIdx));
2167  using SolventModule = BlackOilSolventModule<TypeTag, true>;
2168  if (Fsolgas > SolventModule::cutOff) { // same cutoff as in the solvent model to avoid division by zero
2169  const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
2170  const auto& ssfnKrg = SolventModule::ssfnKrg(satid);
2171  const auto& ssfnKrs = SolventModule::ssfnKrs(satid);
2172  mob[activeGasCompIdx] *= extendEval(ssfnKrg.eval(1-Fsolgas, /*extrapolate=*/true));
2173  mob[Indices::contiSolventEqIdx] = extendEval(ssfnKrs.eval(Fsolgas, /*extrapolate=*/true) * relativePerms[activeGasCompIdx] / intQuants.solventViscosity());
2174  }
2175  }
2176  }
2177 
2178  if (this->isInjector() && !this->inj_fc_multiplier_.empty()) {
2179  const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
2180  const auto& connections = this->well_ecl_.getConnections();
2181  const auto& connection = connections[perf_ecl_index];
2182  if (connection.filterCakeActive()) {
2183  std::ranges::transform(mob, mob.begin(),
2184  [mult = this->inj_fc_multiplier_[local_perf_index]]
2185  (const auto val)
2186  { return val * mult; });
2187  }
2188  }
2189  }
2190 
2191 
2192  template<typename TypeTag>
2193  bool
2194  WellInterface<TypeTag>::
2195  updateWellStateWithTHPTargetProd(const Simulator& simulator,
2196  WellStateType& well_state,
2197  const GroupStateHelperType& groupStateHelper) const
2198  {
2199  auto& deferred_logger = groupStateHelper.deferredLogger();
2200  OPM_TIMEFUNCTION();
2201  const auto& summary_state = simulator.vanguard().summaryState();
2202 
2203  auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
2204  simulator, groupStateHelper, summary_state, this->getALQ(well_state), /*iterate_if_no_solution */ false);
2205  if (bhp_at_thp_limit) {
2206  std::vector<Scalar> rates(this->number_of_phases_, 0.0);
2207  if (thp_update_iterations) {
2208  computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
2209  groupStateHelper, rates);
2210  } else {
2211  computeWellRatesWithBhp(simulator, *bhp_at_thp_limit,
2212  rates, deferred_logger);
2213  }
2214  auto& ws = well_state.well(this->name());
2215  ws.surface_rates = rates;
2216  ws.bhp = *bhp_at_thp_limit;
2217  ws.thp = this->getTHPConstraint(summary_state);
2218  return true;
2219  } else {
2220  return false;
2221  }
2222  }
2223 
2224  template<typename TypeTag>
2225  std::optional<typename WellInterface<TypeTag>::Scalar>
2226  WellInterface<TypeTag>::
2227  computeBhpAtThpLimitProdWithAlqUsingIPR(const Simulator& simulator,
2228  const WellStateType& well_state,
2229  Scalar bhp,
2230  const SummaryState& summary_state,
2231  const Scalar alq_value)
2232  {
2233  OPM_TIMEFUNCTION();
2234  WellStateType well_state_copy = well_state;
2235  const auto& groupStateHelper = simulator.problem().wellModel().groupStateHelper();
2236  GroupStateHelperType groupStateHelper_copy = groupStateHelper;
2237  auto well_guard = groupStateHelper_copy.pushWellState(well_state_copy);
2238  const double dt = simulator.timeStepSize();
2239  const bool converged = this->solveWellWithBhp(
2240  simulator, dt, bhp, groupStateHelper_copy, well_state_copy
2241  );
2242 
2243  bool zero_rates;
2244  auto rates = well_state_copy.well(this->index_of_well_).surface_rates;
2245  zero_rates = true;
2246  for (std::size_t p = 0; p < rates.size(); ++p) {
2247  zero_rates &= rates[p] == 0.0;
2248  }
2249  // For zero rates or unconverged bhp the implicit IPR is problematic.
2250  // Use the old approach for now
2251  if (zero_rates || !converged) {
2252  return this->computeBhpAtThpLimitProdWithAlq(simulator, groupStateHelper_copy, summary_state, alq_value, /*iterate_if_no_solution */ false);
2253  }
2254  this->updateIPRImplicit(simulator, groupStateHelper_copy, well_state_copy);
2255  this->adaptRatesForVFP(rates);
2256  return WellBhpThpCalculator(*this).estimateStableBhp(well_state_copy, this->well_ecl_, rates, this->getRefDensity(), summary_state, alq_value);
2257  }
2258 
2259  template <typename TypeTag>
2260  void
2261  WellInterface<TypeTag>::
2262  computeConnLevelProdInd(const FluidState& fs,
2263  const std::function<Scalar(const Scalar)>& connPICalc,
2264  const std::vector<Scalar>& mobility,
2265  Scalar* connPI) const
2266  {
2267  const int np = this->number_of_phases_;
2268  for (int p = 0; p < np; ++p) {
2269  // Note: E100's notion of PI value phase mobility includes
2270  // the reciprocal FVF.
2271  const int canonical_phase_idx = FluidSystem::activeToCanonicalPhaseIdx(p);
2272  const auto connMob =
2273  mobility[FluidSystem::activePhaseToActiveCompIdx(p)] * fs.invB(canonical_phase_idx).value();
2274 
2275  connPI[p] = connPICalc(connMob);
2276  }
2277 
2278  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2279  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2280  {
2281  const auto io = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2282  const auto ig = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2283 
2284  const auto vapoil = connPI[ig] * fs.Rv().value();
2285  const auto disgas = connPI[io] * fs.Rs().value();
2286 
2287  connPI[io] += vapoil;
2288  connPI[ig] += disgas;
2289  }
2290  }
2291 
2292 
2293  template <typename TypeTag>
2294  void
2295  WellInterface<TypeTag>::
2296  computeConnLevelInjInd(const FluidState& fs,
2297  const Phase preferred_phase,
2298  const std::function<Scalar(const Scalar)>& connIICalc,
2299  const std::vector<Scalar>& mobility,
2300  Scalar* connII,
2301  DeferredLogger& deferred_logger) const
2302  {
2303  auto phase_pos = 0;
2304  if (preferred_phase == Phase::GAS) {
2305  phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
2306  }
2307  else if (preferred_phase == Phase::OIL) {
2308  phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
2309  }
2310  else if (preferred_phase == Phase::WATER) {
2311  phase_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
2312  }
2313  else {
2314  OPM_DEFLOG_THROW(NotImplemented,
2315  fmt::format("Unsupported Injector Type ({}) "
2316  "for well {} during connection I.I. calculation",
2317  static_cast<int>(preferred_phase), this->name()),
2318  deferred_logger);
2319  }
2320 
2321  const auto mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2322  const int canonicalPhaseIdx = FluidSystem::activeToCanonicalPhaseIdx(phase_pos);
2323  connII[phase_pos] = connIICalc(mt * fs.invB(canonicalPhaseIdx).value());
2324  }
2325 
2326  template<typename TypeTag>
2327  template<class GasLiftSingleWell>
2328  std::unique_ptr<GasLiftSingleWell>
2329  WellInterface<TypeTag>::
2330  initializeGliftWellTest_(const Simulator& simulator,
2331  WellStateType& well_state,
2332  const GroupState<Scalar>& group_state,
2333  GLiftEclWells& ecl_well_map,
2334  DeferredLogger& deferred_logger)
2335  {
2336  // Instantiate group info object (without initialization) since it is needed in GasLiftSingleWell
2337  auto& comm = simulator.vanguard().grid().comm();
2338  ecl_well_map.try_emplace(this->name(), &(this->wellEcl()), this->indexOfWell());
2339  const auto& iterCtx = simulator.problem().iterationContext();
2340  GasLiftGroupInfo<Scalar, IndexTraits> group_info {
2341  ecl_well_map,
2342  simulator.vanguard().schedule(),
2343  simulator.vanguard().summaryState(),
2344  simulator.episodeIndex(),
2345  iterCtx,
2346  deferred_logger,
2347  well_state,
2348  group_state,
2349  comm,
2350  false
2351  };
2352 
2353  // Return GasLiftSingleWell object to use the wellTestALQ() function
2354  std::set<int> sync_groups;
2355  const auto& summary_state = simulator.vanguard().summaryState();
2356  return std::make_unique<GasLiftSingleWell>(*this,
2357  simulator,
2358  summary_state,
2359  deferred_logger,
2360  well_state,
2361  group_state,
2362  group_info,
2363  sync_groups,
2364  comm,
2365  false);
2366 
2367  }
2368 
2369 } // namespace Opm
2370 
2371 #endif
Opm::WellInterface::WellInterface
WellInterface(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters &param, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_conservation_quantities, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar >> &perf_data)
Constructor.
Definition: WellInterface_impl.hpp:59
Opm
Structs needed for tpfalinearizer and its gpuparams struct extracted to be defined in one place that ...
Definition: blackoilbioeffectsmodules.hh:45
Opm::WellInterface
Definition: BlackoilWellModelGasLift.hpp:38
Opm::ParallelWellInfo
Class encapsulating some information about parallel wells.
Definition: MSWellHelpers.hpp:34
Opm::WellInterfaceIndices
Definition: MultisegmentWellAssemble.hpp:38
Opm::DeferredLogger
Definition: DeferredLogger.hpp:56
Opm::WellInterface::initializeProducerWellState
void initializeProducerWellState(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) const
Modify the well_state&#39;s rates if there is only one nonzero rate.
Definition: WellInterface_impl.hpp:1822
Opm::PerforationData
Static data associated with a well perforation.
Definition: BlackoilWellModelRestart.hpp:41
Opm::WellState
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: TemperatureModel.hpp:65