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