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