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