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