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