My Project
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/input/eclipse/Schedule/ScheduleTypes.hpp>
23#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24#include <opm/simulators/wells/GroupState.hpp>
25#include <opm/simulators/wells/TargetCalculator.hpp>
26
27#include <dune/common/version.hh>
28
29namespace Opm
30{
31
32
33 template<typename TypeTag>
35 WellInterface(const Well& well,
36 const ParallelWellInfo& pw_info,
37 const int time_step,
38 const ModelParameters& param,
39 const RateConverterType& rate_converter,
40 const int pvtRegionIdx,
41 const int num_components,
42 const int num_phases,
43 const int index_of_well,
44 const std::vector<PerforationData>& perf_data)
45 : WellInterfaceIndices<FluidSystem,Indices,Scalar>(well,
46 pw_info,
47 time_step,
48 rate_converter,
49 pvtRegionIdx,
50 num_components,
51 num_phases,
52 index_of_well,
53 perf_data)
54 , param_(param)
55 {
56 connectionRates_.resize(this->number_of_perforations_);
57
58 if constexpr (has_solvent || has_zFraction) {
59 if (well.isInjector()) {
60 auto injectorType = this->well_ecl_.injectorType();
61 if (injectorType == InjectorType::GAS) {
62 this->wsolvent_ = this->well_ecl_.getSolventFraction();
63 }
64 }
65 }
66 }
67
68
69 template<typename TypeTag>
70 void
72 init(const PhaseUsage* phase_usage_arg,
73 const std::vector<double>& /* depth_arg */,
74 const double gravity_arg,
75 const int /* num_cells */,
76 const std::vector< Scalar >& B_avg,
77 const bool changed_to_open_this_step)
78 {
79 this->phase_usage_ = phase_usage_arg;
80 this->gravity_ = gravity_arg;
81 B_avg_ = B_avg;
82 this->changed_to_open_this_step_ = changed_to_open_this_step;
83 }
84
85
86
87
88 template<typename TypeTag>
89 double
90 WellInterface<TypeTag>::
91 wpolymer() const
92 {
93 if constexpr (has_polymer) {
94 auto injectorType = this->well_ecl_.injectorType();
95
96 if (injectorType == InjectorType::WATER) {
97 WellPolymerProperties polymer = this->well_ecl_.getPolymerProperties();
98 const double polymer_injection_concentration = polymer.m_polymerConcentration;
99 return polymer_injection_concentration;
100 } else {
101 // Not a water injection well => no polymer.
102 return 0.0;
103 }
104 }
105
106 return 0.0;
107 }
108
109
110
111
112
113 template<typename TypeTag>
114 double
115 WellInterface<TypeTag>::
116 wfoam() const
117 {
118 if constexpr (has_foam) {
119 auto injectorType = this->well_ecl_.injectorType();
120
121 if (injectorType == InjectorType::GAS) {
122 WellFoamProperties fprop = this->well_ecl_.getFoamProperties();
123 return fprop.m_foamConcentration;
124 } else {
125 // Not a gas injection well => no foam.
126 return 0.0;
127 }
128 }
129
130 return 0.0;
131 }
132
133
134
135 template<typename TypeTag>
136 double
137 WellInterface<TypeTag>::
138 wsalt() const
139 {
140 if constexpr (has_brine) {
141 auto injectorType = this->well_ecl_.injectorType();
142
143 if (injectorType == InjectorType::WATER) {
144 WellBrineProperties fprop = this->well_ecl_.getBrineProperties();
145 return fprop.m_saltConcentration;
146 } else {
147 // Not a water injection well => no salt (?).
148 return 0.0;
149 }
150 }
151
152 return 0.0;
153 }
154
155 template<typename TypeTag>
156 double
157 WellInterface<TypeTag>::
158 wmicrobes() const
159 {
160 if constexpr (has_micp) {
161 auto injectorType = this->well_ecl_.injectorType();
162
163 if (injectorType == InjectorType::WATER) {
164 WellMICPProperties microbes = this->well_ecl_.getMICPProperties();
165 const double microbial_injection_concentration = microbes.m_microbialConcentration;
166 return microbial_injection_concentration;
167 } else {
168 // Not a water injection well => no microbes.
169 return 0.0;
170 }
171 }
172
173 return 0.0;
174 }
175
176 template<typename TypeTag>
177 double
178 WellInterface<TypeTag>::
179 woxygen() const
180 {
181 if constexpr (has_micp) {
182 auto injectorType = this->well_ecl_.injectorType();
183
184 if (injectorType == InjectorType::WATER) {
185 WellMICPProperties oxygen = this->well_ecl_.getMICPProperties();
186 const double oxygen_injection_concentration = oxygen.m_oxygenConcentration;
187 return oxygen_injection_concentration;
188 } else {
189 // Not a water injection well => no oxygen.
190 return 0.0;
191 }
192 }
193
194 return 0.0;
195 }
196
197 // The urea injection concentration is scaled down by a factor of 10, since its value
198 // can be much bigger than 1 (not doing this slows the simulations). The
199 // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing
200 // the reactions and also when writing the output files (vtk and eclipse format, i.e.,
201 // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively).
202
203 template<typename TypeTag>
204 double
205 WellInterface<TypeTag>::
206 wurea() const
207 {
208 if constexpr (has_micp) {
209 auto injectorType = this->well_ecl_.injectorType();
210
211 if (injectorType == InjectorType::WATER) {
212 WellMICPProperties urea = this->well_ecl_.getMICPProperties();
213 const double urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
214 return urea_injection_concentration;
215 } else {
216 // Not a water injection well => no urea.
217 return 0.0;
218 }
219 }
220
221 return 0.0;
222 }
223
224 template<typename TypeTag>
225 bool
226 WellInterface<TypeTag>::
227 updateWellControl(const Simulator& ebos_simulator,
228 const IndividualOrGroup iog,
229 WellState& well_state,
230 const GroupState& group_state,
231 DeferredLogger& deferred_logger) /* const */
232 {
233 if (this->wellIsStopped()) {
234 return false;
235 }
236
237 const auto& summaryState = ebos_simulator.vanguard().summaryState();
238 const auto& schedule = ebos_simulator.vanguard().schedule();
239 const auto& well = this->well_ecl_;
240 auto& ws = well_state.well(this->index_of_well_);
241 std::string from;
242 if (well.isInjector()) {
243 from = Well::InjectorCMode2String(ws.injection_cmode);
244 } else {
245 from = Well::ProducerCMode2String(ws.production_cmode);
246 }
247 bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
248
249 if (oscillating) {
250 // only output frist time
251 bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
252 if (output) {
253 std::ostringstream ss;
254 ss << " The control model for well " << this->name()
255 << " is oscillating\n"
256 << " We don't allow for more than "
257 << param_.max_number_of_well_switches_
258 << " switches. The control is kept at " << from;
259 deferred_logger.info(ss.str());
260 // add one more to avoid outputting the same info again
261 this->well_control_log_.push_back(from);
262 }
263 return false;
264 }
265 bool changed = false;
266 if (iog == IndividualOrGroup::Individual) {
267 changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger);
268 } else if (iog == IndividualOrGroup::Group) {
269 changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
270 } else {
271 assert(iog == IndividualOrGroup::Both);
272 changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
273 }
274 Parallel::Communication cc = ebos_simulator.vanguard().grid().comm();
275 // checking whether control changed
276 if (changed) {
277 std::string to;
278 if (well.isInjector()) {
279 to = Well::InjectorCMode2String(ws.injection_cmode);
280 } else {
281 to = Well::ProducerCMode2String(ws.production_cmode);
282 }
283 std::ostringstream ss;
284 ss << " Switching control mode for well " << this->name()
285 << " from " << from
286 << " to " << to;
287 if (cc.size() > 1) {
288 ss << " on rank " << cc.rank();
289 }
290 deferred_logger.debug(ss.str());
291
292 this->well_control_log_.push_back(from);
293 updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
294 updatePrimaryVariables(well_state, deferred_logger);
295 }
296
297 return changed;
298 }
299
300
301
302 template<typename TypeTag>
303 void
304 WellInterface<TypeTag>::
305 wellTesting(const Simulator& simulator,
306 const double simulation_time,
307 /* const */ WellState& well_state,
308 const GroupState& group_state,
309 WellTestState& well_test_state,
310 DeferredLogger& deferred_logger)
311 {
312 deferred_logger.info(" well " + this->name() + " is being tested");
313
314 WellState well_state_copy = well_state;
315 auto& ws = well_state_copy.well(this->indexOfWell());
316
317 updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
318 calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
319 updatePrimaryVariables(well_state_copy, deferred_logger);
320 initPrimaryVariablesEvaluation();
321
322 if (this->isProducer()) {
323 gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
324 }
325
326 WellTestState welltest_state_temp;
327
328 bool testWell = true;
329 // if a well is closed because all completions are closed, we need to check each completion
330 // individually. We first open all completions, then we close one by one by calling updateWellTestState
331 // untill the number of closed completions do not increase anymore.
332 while (testWell) {
333 const size_t original_number_closed_completions = welltest_state_temp.num_closed_completions();
334 bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger);
335 if (!converged) {
336 const auto msg = fmt::format("WTEST: Well {} is not solvable (physical)", this->name());
337 deferred_logger.debug(msg);
338 return;
339 }
340
341
342 updateWellOperability(simulator, well_state_copy, deferred_logger);
343 if ( !this->isOperableAndSolvable() ) {
344 const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
345 deferred_logger.debug(msg);
346 return;
347 }
348
349 std::vector<double> potentials;
350 try {
351 computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
352 } catch (const std::exception& e) {
353 const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during testing for re-opening: ") + e.what();
354 deferred_logger.info(msg);
355 return;
356 }
357 const int np = well_state_copy.numPhases();
358 for (int p = 0; p < np; ++p) {
359 ws.well_potentials[p] = std::max(0.0, potentials[p]);
360 }
361 this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
362 this->closeCompletions(welltest_state_temp);
363
364 // Stop testing if the well is closed or shut due to all completions shut
365 // Also check if number of completions has increased. If the number of closed completions do not increased
366 // we stop the testing.
367 // TODO: it can be tricky here, if the well is shut/closed due to other reasons
368 if ( welltest_state_temp.num_closed_wells() > 0 ||
369 (original_number_closed_completions == welltest_state_temp.num_closed_completions()) ) {
370 testWell = false; // this terminates the while loop
371 }
372 }
373
374 // update wellTestState if the well test succeeds
375 if (!welltest_state_temp.well_is_closed(this->name())) {
376 well_test_state.open_well(this->name());
377
378 std::string msg = std::string("well ") + this->name() + std::string(" is re-opened");
379 deferred_logger.info(msg);
380
381 // also reopen completions
382 for (auto& completion : this->well_ecl_.getCompletions()) {
383 if (!welltest_state_temp.completion_is_closed(this->name(), completion.first))
384 well_test_state.open_completion(this->name(), completion.first);
385 }
386 // set the status of the well_state to open
387 ws.open();
388 well_state = well_state_copy;
389 }
390 }
391
392
393
394
395 template<typename TypeTag>
396 bool
397 WellInterface<TypeTag>::
398 iterateWellEquations(const Simulator& ebosSimulator,
399 const double dt,
400 WellState& well_state,
401 const GroupState& group_state,
402 DeferredLogger& deferred_logger)
403 {
404 const auto& summary_state = ebosSimulator.vanguard().summaryState();
405 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
406 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
407 bool converged = false;
408 try {
409 converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
410 } catch (NumericalIssue& e ) {
411 const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. ";
412 deferred_logger.warning("INNER_ITERATION_FAILED", msg);
413 converged = false;
414 }
415 return converged;
416 }
417
418
419 template<typename TypeTag>
420 bool
421 WellInterface<TypeTag>::
422 solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
423 DeferredLogger& deferred_logger)
424 {
425 // keep a copy of the original well state
426 const WellState well_state0 = well_state;
427 const double dt = ebosSimulator.timeStepSize();
428 const auto& summary_state = ebosSimulator.vanguard().summaryState();
429 const bool has_thp_limit = this->wellHasTHPConstraints(summary_state);
430 if (has_thp_limit)
431 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
432 else
433 well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
434
435 const bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
436 if (converged) {
437 deferred_logger.debug("WellTest: Well equation for well " + this->name() + " converged");
438 return true;
439 }
440 const int max_iter = param_.max_welleq_iter_;
441 deferred_logger.debug("WellTest: Well equation for well " + this->name() + " failed converging in "
442 + std::to_string(max_iter) + " iterations");
443 well_state = well_state0;
444 return false;
445 }
446
447
448 template<typename TypeTag>
449 void
450 WellInterface<TypeTag>::
451 solveWellEquation(const Simulator& ebosSimulator,
452 WellState& well_state,
453 const GroupState& group_state,
454 DeferredLogger& deferred_logger)
455 {
456 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
457 return;
458
459 // keep a copy of the original well state
460 const WellState well_state0 = well_state;
461 const double dt = ebosSimulator.timeStepSize();
462 bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
463
464 // Newly opened wells with THP control sometimes struggles to
465 // converge due to bad initial guess. Or due to the simple fact
466 // that the well needs to change to another control.
467 // We therefore try to solve the well with BHP control to get
468 // an better initial guess.
469 // If the well is supposed to operate under THP control
470 // "updateWellControl" will switch it back to THP later.
471 if (!converged) {
472 auto& ws = well_state.well(this->indexOfWell());
473 bool thp_control = false;
474 if (this->well_ecl_.isInjector()) {
475 thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
476 if (thp_control) {
477 ws.injection_cmode = Well::InjectorCMode::BHP;
478 this->well_control_log_.push_back(Well::InjectorCMode2String(Well::InjectorCMode::THP));
479 }
480 } else {
481 thp_control = ws.production_cmode == Well::ProducerCMode::THP;
482 if (thp_control) {
483 ws.production_cmode = Well::ProducerCMode::BHP;
484 this->well_control_log_.push_back(Well::ProducerCMode2String(Well::ProducerCMode::THP));
485 }
486 }
487 if (thp_control) {
488 const std::string msg = std::string("The newly opened well ") + this->name()
489 + std::string(" with THP control did not converge during inner iterations, we try again with bhp control");
490 deferred_logger.debug(msg);
491 converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
492 }
493 }
494
495 if (!converged) {
496 const int max_iter = param_.max_welleq_iter_;
497 deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
498 + std::to_string(max_iter) + " iterations");
499 well_state = well_state0;
500 }
501 }
502
503
504
505 template <typename TypeTag>
506 void
507 WellInterface<TypeTag>::
508 assembleWellEq(const Simulator& ebosSimulator,
509 const double dt,
510 WellState& well_state,
511 const GroupState& group_state,
512 DeferredLogger& deferred_logger)
513 {
514 const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
515
516 if (param_.check_well_operability_iter_)
517 checkWellOperability(ebosSimulator, well_state, deferred_logger);
518
519 // only use inner well iterations for the first newton iterations.
520 const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
521 if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
522 this->operability_status_.solvable = true;
523 bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
524
525 // unsolvable wells are treated as not operable and will not be solved for in this iteration.
526 if (!converged) {
527 if (param_.shut_unsolvable_wells_)
528 this->operability_status_.solvable = false;
529 }
530 }
531 if (this->operability_status_.has_negative_potentials) {
532 auto well_state_copy = well_state;
533 std::vector<double> potentials;
534 try {
535 computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
536 } catch (const std::exception& e) {
537 const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
538 deferred_logger.info(msg);
539 this->operability_status_.has_negative_potentials = true;
540 }
541 auto& ws = well_state.well(this->indexOfWell());
542 const int np = well_state.numPhases();
543 for (int p = 0; p < np; ++p) {
544 ws.well_potentials[p] = std::max(0.0, potentials[p]);
545 }
546 }
547 this->changed_to_open_this_step_ = false;
548 const bool well_operable = this->operability_status_.isOperableAndSolvable();
549
550 if (!well_operable && old_well_operable) {
551 if (this->well_ecl_.getAutomaticShutIn()) {
552 deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
553 } else {
554 if (!this->wellIsStopped()) {
555 deferred_logger.info(" well " + this->name() + " gets STOPPED during iteration ");
556 this->stopWell();
557 changed_to_stopped_this_step_ = true;
558 }
559 }
560 } else if (well_operable && !old_well_operable) {
561 deferred_logger.info(" well " + this->name() + " gets REVIVED during iteration ");
562 this->openWell();
563 changed_to_stopped_this_step_ = false;
564 this->changed_to_open_this_step_ = true;
565 }
566
567 const auto& summary_state = ebosSimulator.vanguard().summaryState();
568 const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
569 const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
570 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
571 }
572
573 template<typename TypeTag>
574 bool
575 WellInterface<TypeTag>::isPressureControlled(const WellState& well_state) const
576 {
577 bool thp_controlled_well = false;
578 bool bhp_controlled_well = false;
579 const auto& ws = well_state.well(this->index_of_well_);
580 if (this->isInjector()) {
581 const Well::InjectorCMode& current = ws.injection_cmode;
582 if (current == Well::InjectorCMode::THP) {
583 thp_controlled_well = true;
584 }
585 if (current == Well::InjectorCMode::BHP) {
586 bhp_controlled_well = true;
587 }
588 } else {
589 const Well::ProducerCMode& current = ws.production_cmode;
590 if (current == Well::ProducerCMode::THP) {
591 thp_controlled_well = true;
592 }
593 if (current == Well::ProducerCMode::BHP) {
594 bhp_controlled_well = true;
595 }
596 }
597 bool ispressureControlled = (bhp_controlled_well || thp_controlled_well);
598 return ispressureControlled;
599 }
600
601 template<typename TypeTag>
602 void
603 WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const
604 {
605 if(!this->isOperableAndSolvable() && !this->wellIsStopped())
606 return;
607
608 for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
609 if (this->cells()[perfIdx] == cellIdx) {
610 for (int i = 0; i < RateVector::dimension; ++i) {
611 rates[i] += connectionRates_[perfIdx][i];
612 }
613 }
614 }
615 }
616
617 template<typename TypeTag>
618 typename WellInterface<TypeTag>::Scalar
619 WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const {
620 for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
621 if (this->cells()[perfIdx] == cellIdx) {
622 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
623 return connectionRates_[perfIdx][activeCompIdx].value();
624 }
625 }
626 // this is not thread safe
627 OPM_THROW(std::invalid_argument, "The well with name " + this->name()
628 + " does not perforate cell " + std::to_string(cellIdx));
629 return 0.0;
630 }
631
632
633
634
635 template<typename TypeTag>
636 void
637 WellInterface<TypeTag>::
638 checkWellOperability(const Simulator& ebos_simulator,
639 const WellState& well_state,
640 DeferredLogger& deferred_logger)
641 {
642
643 if (!param_.check_well_operability_) {
644 return;
645 }
646
647 if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
648 return;
649 }
650
651 updateWellOperability(ebos_simulator, well_state, deferred_logger);
652 if (!this->operability_status_.isOperableAndSolvable()) {
653 this->operability_status_.use_vfpexplicit = true;
654 deferred_logger.debug("EXPLICIT_LOOKUP_VFP",
655 "well not operable, trying with explicit vfp lookup: " + this->name());
656 updateWellOperability(ebos_simulator, well_state, deferred_logger);
657 }
658 }
659
660 template<typename TypeTag>
661 void
662 WellInterface<TypeTag>::
663 gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
664 WellState& well_state,
665 DeferredLogger& deferred_logger)
666 {
667 const auto& summary_state = ebos_simulator.vanguard().summaryState();
668 const auto& well_name = this->name();
669 if (!this->wellHasTHPConstraints(summary_state)) {
670 const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
671 deferred_logger.info(msg);
672 return;
673 }
674 const auto& well_ecl = this->wellEcl();
675 const auto& schedule = ebos_simulator.vanguard().schedule();
676 auto report_step_idx = ebos_simulator.episodeIndex();
677 const auto& glo = schedule.glo(report_step_idx);
678 if (!glo.has_well(well_name)) {
679 const std::string msg = fmt::format(
680 "GLIFT WTEST: Well {} : Gas Lift not activated: "
681 "WLIFTOPT is probably missing. Skipping.", well_name);
682 deferred_logger.info(msg);
683 return;
684 }
685 const auto& gl_well = glo.well(well_name);
686 auto& max_alq_optional = gl_well.max_rate();
687 double max_alq;
688 if (max_alq_optional) {
689 max_alq = *max_alq_optional;
690 }
691 else {
692 const auto& controls = well_ecl.productionControls(summary_state);
693 const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
694 const auto& alq_values = table.getALQAxis();
695 max_alq = alq_values.back();
696 }
697 well_state.setALQ(well_name, max_alq);
698 const std::string msg = fmt::format(
699 "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
700 well_name, max_alq);
701 deferred_logger.info(msg);
702 }
703
704 template<typename TypeTag>
705 void
706 WellInterface<TypeTag>::
707 updateWellOperability(const Simulator& ebos_simulator,
708 const WellState& well_state,
709 DeferredLogger& deferred_logger)
710 {
711 this->operability_status_.resetOperability();
712
713 bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
714 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
715 bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
716 well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
717
718 // Operability checking is not free
719 // Only check wells under BHP and THP control
720 bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
721 if (check_thp || bhp_controlled) {
722 updateIPR(ebos_simulator, deferred_logger);
723 checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
724 }
725 // we do some extra checking for wells under THP control.
726 if (check_thp) {
727 checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
728 }
729 }
730
731
732 template<typename TypeTag>
733 void
734 WellInterface<TypeTag>::
735 updateWellStateWithTarget(const Simulator& ebos_simulator,
736 const GroupState& group_state,
737 WellState& well_state,
738 DeferredLogger& deferred_logger) const
739 {
740
741 // only bhp and wellRates are used to initilize the primaryvariables for standard wells
742 const auto& well = this->well_ecl_;
743 const int well_index = this->index_of_well_;
744 auto& ws = well_state.well(well_index);
745 const auto& pu = this->phaseUsage();
746 const int np = well_state.numPhases();
747 const auto& summaryState = ebos_simulator.vanguard().summaryState();
748 const auto& schedule = ebos_simulator.vanguard().schedule();
749
750 if (this->wellIsStopped()) {
751 for (int p = 0; p<np; ++p) {
752 ws.surface_rates[p] = 0;
753 }
754 ws.thp = 0;
755 return;
756 }
757
758 if (this->isInjector() )
759 {
760 const auto& controls = well.injectionControls(summaryState);
761
762 InjectorType injectorType = controls.injector_type;
763 int phasePos;
764 switch (injectorType) {
765 case InjectorType::WATER:
766 {
767 phasePos = pu.phase_pos[BlackoilPhases::Aqua];
768 break;
769 }
770 case InjectorType::OIL:
771 {
772 phasePos = pu.phase_pos[BlackoilPhases::Liquid];
773 break;
774 }
775 case InjectorType::GAS:
776 {
777 phasePos = pu.phase_pos[BlackoilPhases::Vapour];
778 break;
779 }
780 default:
781 OPM_DEFLOG_THROW(std::runtime_error, "Expected WATER, OIL or GAS as type for injectors " + this->name(), deferred_logger );
782 }
783
784 const auto current = ws.injection_cmode;
785
786 switch(current) {
787 case Well::InjectorCMode::RATE:
788 {
789 ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
790 if(this->rsRvInj() > 0) {
791 if (injectorType == InjectorType::OIL && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
792 ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.surface_rate * this->rsRvInj();
793 } else if (injectorType == InjectorType::GAS && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
794 ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.surface_rate * this->rsRvInj();
795 } else {
796 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 );
797 }
798 }
799 break;
800 }
801
802 case Well::InjectorCMode::RESV:
803 {
804 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
805 this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
806 const double coeff = convert_coeff[phasePos];
807 ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
808 break;
809 }
810
811 case Well::InjectorCMode::THP:
812 {
813 auto rates = ws.surface_rates;
814 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
815 ws.bhp = bhp;
816 ws.thp = this->getTHPConstraint(summaryState);
817
818 // if the total rates are negative or zero
819 // we try to provide a better intial well rate
820 // using the well potentials
821 double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
822 if (total_rate <= 0.0)
823 ws.surface_rates = ws.well_potentials;
824
825 break;
826 }
827 case Well::InjectorCMode::BHP:
828 {
829 ws.bhp = controls.bhp_limit;
830 double total_rate = 0.0;
831 for (int p = 0; p<np; ++p) {
832 total_rate += ws.surface_rates[p];
833 }
834 // if the total rates are negative or zero
835 // we try to provide a better intial well rate
836 // using the well potentials
837 if (total_rate <= 0.0)
838 ws.surface_rates = ws.well_potentials;
839
840 break;
841 }
842 case Well::InjectorCMode::GRUP:
843 {
844 assert(well.isAvailableForGroupControl());
845 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
846 const double efficiencyFactor = well.getEfficiencyFactor();
847 std::optional<double> target =
848 this->getGroupInjectionTargetRate(group,
849 well_state,
850 group_state,
851 schedule,
852 summaryState,
853 injectorType,
854 efficiencyFactor,
855 deferred_logger);
856 if (target)
857 ws.surface_rates[phasePos] = *target;
858 break;
859 }
860 case Well::InjectorCMode::CMODE_UNDEFINED:
861 {
862 OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger );
863 }
864
865 }
866 // for wells with zero injection rate, if we assign exactly zero rate,
867 // we will have to assume some trivial composition in the wellbore.
868 // here, we use some small value (about 0.01 m^3/day ~= 1.e-7) to initialize
869 // the zero rate target, then we can use to retain the composition information
870 // within the wellbore from the previous result, and hopefully it is a good
871 // initial guess for the zero rate target.
872 ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
873 }
874 //Producer
875 else
876 {
877 const auto current = ws.production_cmode;
878 const auto& controls = well.productionControls(summaryState);
879 switch (current) {
880 case Well::ProducerCMode::ORAT:
881 {
882 double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
883 // for trivial rates or opposite direction we don't just scale the rates
884 // but use either the potentials or the mobility ratio to initial the well rates
885 if (current_rate > 0.0) {
886 for (int p = 0; p<np; ++p) {
887 ws.surface_rates[p] *= controls.oil_rate/current_rate;
888 }
889 } else {
890 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
891 double control_fraction = fractions[pu.phase_pos[Oil]];
892 if (control_fraction != 0.0) {
893 for (int p = 0; p<np; ++p) {
894 ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
895 }
896 }
897 }
898 break;
899 }
900 case Well::ProducerCMode::WRAT:
901 {
902 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
903 // for trivial rates or opposite direction we don't just scale the rates
904 // but use either the potentials or the mobility ratio to initial the well rates
905 if (current_rate > 0.0) {
906 for (int p = 0; p<np; ++p) {
907 ws.surface_rates[p] *= controls.water_rate/current_rate;
908 }
909 } else {
910 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
911 double control_fraction = fractions[pu.phase_pos[Water]];
912 if (control_fraction != 0.0) {
913 for (int p = 0; p<np; ++p) {
914 ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
915 }
916 }
917 }
918 break;
919 }
920 case Well::ProducerCMode::GRAT:
921 {
922 double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
923 // or trivial rates or opposite direction we don't just scale the rates
924 // but use either the potentials or the mobility ratio to initial the well rates
925 if (current_rate > 0.0) {
926 for (int p = 0; p<np; ++p) {
927 ws.surface_rates[p] *= controls.gas_rate/current_rate;
928 }
929 } else {
930 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
931 double control_fraction = fractions[pu.phase_pos[Gas]];
932 if (control_fraction != 0.0) {
933 for (int p = 0; p<np; ++p) {
934 ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
935 }
936 }
937 }
938
939 break;
940
941 }
942 case Well::ProducerCMode::LRAT:
943 {
944 double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
945 - ws.surface_rates[ pu.phase_pos[Oil] ];
946 // or trivial rates or opposite direction we don't just scale the rates
947 // but use either the potentials or the mobility ratio to initial the well rates
948 if (current_rate > 0.0) {
949 for (int p = 0; p<np; ++p) {
950 ws.surface_rates[p] *= controls.liquid_rate/current_rate;
951 }
952 } else {
953 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
954 double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
955 if (control_fraction != 0.0) {
956 for (int p = 0; p<np; ++p) {
957 ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
958 }
959 }
960 }
961 break;
962 }
963 case Well::ProducerCMode::CRAT:
964 {
965 OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << this->name(), deferred_logger);
966 }
967 case Well::ProducerCMode::RESV:
968 {
969 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
970 this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
971 double total_res_rate = 0.0;
972 for (int p = 0; p<np; ++p) {
973 total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
974 }
975 if (controls.prediction_mode) {
976 // or trivial rates or opposite direction we don't just scale the rates
977 // but use either the potentials or the mobility ratio to initial the well rates
978 if (total_res_rate > 0.0) {
979 for (int p = 0; p<np; ++p) {
980 ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
981 }
982 } else {
983 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
984 for (int p = 0; p<np; ++p) {
985 ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
986 }
987 }
988 } else {
989 std::vector<double> hrates(this->number_of_phases_,0.);
990 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
991 hrates[pu.phase_pos[Water]] = controls.water_rate;
992 }
993 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
994 hrates[pu.phase_pos[Oil]] = controls.oil_rate;
995 }
996 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
997 hrates[pu.phase_pos[Gas]] = controls.gas_rate;
998 }
999 std::vector<double> hrates_resv(this->number_of_phases_,0.);
1000 this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
1001 double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
1002 // or trivial rates or opposite direction we don't just scale the rates
1003 // but use either the potentials or the mobility ratio to initial the well rates
1004 if (total_res_rate > 0.0) {
1005 for (int p = 0; p<np; ++p) {
1006 ws.surface_rates[p] *= target/total_res_rate;
1007 }
1008 } else {
1009 const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
1010 for (int p = 0; p<np; ++p) {
1011 ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
1012 }
1013 }
1014
1015 }
1016 break;
1017 }
1018 case Well::ProducerCMode::BHP:
1019 {
1020 ws.bhp = controls.bhp_limit;
1021 double total_rate = 0.0;
1022 for (int p = 0; p<np; ++p) {
1023 total_rate -= ws.surface_rates[p];
1024 }
1025 // if the total rates are negative or zero
1026 // we try to provide a better intial well rate
1027 // using the well potentials
1028 if (total_rate <= 0.0){
1029 for (int p = 0; p<np; ++p) {
1030 ws.surface_rates[p] = -ws.well_potentials[p];
1031 }
1032 }
1033 break;
1034 }
1035 case Well::ProducerCMode::THP:
1036 {
1037 auto rates = ws.surface_rates;
1038 this->adaptRatesForVFP(rates);
1039 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
1040 ws.bhp = bhp;
1041 ws.thp = this->getTHPConstraint(summaryState);
1042
1043 // if the total rates are negative or zero
1044 // we try to provide a better intial well rate
1045 // using the well potentials
1046 double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
1047 if (total_rate <= 0.0){
1048 for (int p = 0; p<np; ++p) {
1049 ws.surface_rates[p] = -ws.well_potentials[p];
1050 }
1051 }
1052 break;
1053 }
1054 case Well::ProducerCMode::GRUP:
1055 {
1056 assert(well.isAvailableForGroupControl());
1057 const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
1058 const double efficiencyFactor = well.getEfficiencyFactor();
1059 double scale = this->getGroupProductionTargetRate(group,
1060 well_state,
1061 group_state,
1062 schedule,
1063 summaryState,
1064 efficiencyFactor);
1065
1066 // we don't want to scale with zero and get zero rates.
1067 if (scale > 0) {
1068 for (int p = 0; p<np; ++p) {
1069 ws.surface_rates[p] *= scale;
1070 }
1071 ws.trivial_target = false;
1072 } else {
1073 ws.trivial_target = true;
1074 }
1075 break;
1076 }
1077 case Well::ProducerCMode::CMODE_UNDEFINED:
1078 case Well::ProducerCMode::NONE:
1079 {
1080 OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name() , deferred_logger);
1081 }
1082
1083 break;
1084 } // end of switch
1085 }
1086 }
1087
1088 template<typename TypeTag>
1089 std::vector<double>
1090 WellInterface<TypeTag>::
1091 initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const
1092 {
1093 const int np = this->number_of_phases_;
1094 std::vector<double> scaling_factor(np);
1095 const auto& ws = well_state.well(this->index_of_well_);
1096
1097 double total_potentials = 0.0;
1098 for (int p = 0; p<np; ++p) {
1099 total_potentials += ws.well_potentials[p];
1100 }
1101 if (total_potentials > 0) {
1102 for (int p = 0; p<np; ++p) {
1103 scaling_factor[p] = ws.well_potentials[p] / total_potentials;
1104 }
1105 return scaling_factor;
1106 }
1107 // if we don't have any potentials we weight it using the mobilites
1108 // We only need approximation so we don't bother with the vapporized oil and dissolved gas
1109 double total_tw = 0;
1110 const int nperf = this->number_of_perforations_;
1111 for (int perf = 0; perf < nperf; ++perf) {
1112 total_tw += this->well_index_[perf];
1113 }
1114 for (int perf = 0; perf < nperf; ++perf) {
1115 const int cell_idx = this->well_cells_[perf];
1116 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1117 const auto& fs = intQuants.fluidState();
1118 const double well_tw_fraction = this->well_index_[perf] / total_tw;
1119 double total_mobility = 0.0;
1120 for (int p = 0; p < np; ++p) {
1121 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1122 total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1123 }
1124 for (int p = 0; p < np; ++p) {
1125 int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1126 scaling_factor[p] += well_tw_fraction * fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1127 }
1128 }
1129 return scaling_factor;
1130 }
1131
1132
1133
1134 template <typename TypeTag>
1135 void
1137 updateWellStateRates(const Simulator& ebosSimulator,
1138 WellState& well_state,
1139 DeferredLogger& deferred_logger) const
1140 {
1141 // Check if the rates of this well only are single-phase, do nothing
1142 // if more than one nonzero rate.
1143 auto& ws = well_state.well(this->index_of_well_);
1144 int nonzero_rate_index = -1;
1145 const double floating_point_error_epsilon = 1e-14;
1146 for (int p = 0; p < this->number_of_phases_; ++p) {
1147 if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
1148 if (nonzero_rate_index == -1) {
1149 nonzero_rate_index = p;
1150 } else {
1151 // More than one nonzero rate.
1152 return;
1153 }
1154 }
1155 }
1156
1157 // Calculate the rates that follow from the current primary variables.
1158 std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
1159
1160 if (nonzero_rate_index == -1) {
1161 // No nonzero rates.
1162 // Use the computed rate directly
1163 for (int p = 0; p < this->number_of_phases_; ++p) {
1164 ws.surface_rates[p] = well_q_s[this->flowPhaseToEbosCompIdx(p)];
1165 }
1166 return;
1167 }
1168
1169 // Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
1170 const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
1171 const int comp_idx_nz = this->flowPhaseToEbosCompIdx(nonzero_rate_index);
1172 for (int p = 0; p < this->number_of_phases_; ++p) {
1173 if (p != nonzero_rate_index) {
1174 const int comp_idx = this->flowPhaseToEbosCompIdx(p);
1175 double& rate = ws.surface_rates[p];
1176 rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
1177 }
1178 }
1179 }
1180 template<typename TypeTag>
1181 typename WellInterface<TypeTag>::Eval
1182 WellInterface<TypeTag>::getPerfCellPressure(const typename WellInterface<TypeTag>::FluidState& fs) const
1183 {
1184 Eval pressure;
1185 if (Indices::oilEnabled) {
1186 pressure = fs.pressure(FluidSystem::oilPhaseIdx);
1187 } else {
1188 if (Indices::waterEnabled) {
1189 pressure = fs.pressure(FluidSystem::waterPhaseIdx);
1190 } else {
1191 pressure = fs.pressure(FluidSystem::gasPhaseIdx);
1192 }
1193 }
1194 return pressure;
1195 }
1196} // namespace Opm
Definition: DeferredLogger.hpp:57
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:243
Definition: WellInterfaceIndices.hpp:35
Definition: WellInterface.hpp:72
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:1137
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:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
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:327
Definition: BlackoilPhases.hpp:46