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