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