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