StandardWell_impl.hpp
Go to the documentation of this file.
1/*
2 Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3 Copyright 2017 Statoil ASA.
4 Copyright 2016 - 2017 IRIS AS.
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#ifndef OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
23#define OPM_STANDARDWELL_IMPL_HEADER_INCLUDED
24
25// Improve IDE experience
26#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
27#include <config.h>
29#endif
30
31#include <opm/common/Exceptions.hpp>
32
33#include <opm/input/eclipse/Units/Units.hpp>
34
40
41#include <algorithm>
42#include <cstddef>
43#include <functional>
44
45#include <fmt/format.h>
46
47namespace Opm
48{
49
50 template<typename TypeTag>
52 StandardWell(const Well& well,
53 const ParallelWellInfo<Scalar>& pw_info,
54 const int time_step,
55 const ModelParameters& param,
56 const RateConverterType& rate_converter,
57 const int pvtRegionIdx,
58 const int num_conservation_quantities,
59 const int num_phases,
60 const int index_of_well,
61 const std::vector<PerforationData<Scalar>>& perf_data)
62 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_conservation_quantities, num_phases, index_of_well, perf_data)
63 , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices>&>(*this))
64 , regularize_(false)
65 {
67 }
68
69
70
71
72
73 template<typename TypeTag>
74 void
76 init(const std::vector<Scalar>& depth_arg,
77 const Scalar gravity_arg,
78 const std::vector< Scalar >& B_avg,
79 const bool changed_to_open_this_step)
80 {
81 Base::init(depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
82 this->StdWellEval::init(this->perf_depth_, depth_arg, Base::has_polymermw);
83 }
84
85
86
87
88
89 template<typename TypeTag>
90 template<class Value>
91 void
94 const std::vector<Value>& mob,
95 const Value& bhp,
96 const std::vector<Scalar>& Tw,
97 const int perf,
98 const bool allow_cf,
99 std::vector<Value>& cq_s,
100 PerforationRates<Scalar>& perf_rates,
101 DeferredLogger& deferred_logger) const
102 {
103 auto obtain = [this](const Eval& value)
104 {
105 if constexpr (std::is_same_v<Value, Scalar>) {
106 static_cast<void>(this); // suppress clang warning
107 return getValue(value);
108 } else {
109 return this->extendEval(value);
110 }
111 };
112 auto obtainN = [](const auto& value)
113 {
114 if constexpr (std::is_same_v<Value, Scalar>) {
115 return getValue(value);
116 } else {
117 return value;
118 }
119 };
120 auto zeroElem = [this]()
121 {
122 if constexpr (std::is_same_v<Value, Scalar>) {
123 static_cast<void>(this); // suppress clang warning
124 return 0.0;
125 } else {
126 return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
127 }
128 };
129
130 const auto& fs = intQuants.fluidState();
131 const Value pressure = obtain(this->getPerfCellPressure(fs));
132 const Value rs = obtain(fs.Rs());
133 const Value rv = obtain(fs.Rv());
134 const Value rvw = obtain(fs.Rvw());
135 const Value rsw = obtain(fs.Rsw());
136
137 std::vector<Value> b_perfcells_dense(this->numConservationQuantities(), zeroElem());
138 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
139 if (!FluidSystem::phaseIsActive(phaseIdx)) {
140 continue;
141 }
142 const unsigned compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
143 b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
144 }
145 if constexpr (has_solvent) {
146 b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
147 }
148
149 if constexpr (has_zFraction) {
150 if (this->isInjector()) {
151 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
152 b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
153 b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
154 }
155 }
156
157 Value skin_pressure = zeroElem();
158 if (has_polymermw) {
159 if (this->isInjector()) {
160 const int pskin_index = Bhp + 1 + this->numLocalPerfs() + perf;
161 skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
162 }
163 }
164
165 // surface volume fraction of fluids within wellbore
166 std::vector<Value> cmix_s(this->numConservationQuantities(), zeroElem());
167 for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
168 cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
169 }
170
171 computePerfRate(mob,
172 pressure,
173 bhp,
174 rs,
175 rv,
176 rvw,
177 rsw,
178 b_perfcells_dense,
179 Tw,
180 perf,
181 allow_cf,
182 skin_pressure,
183 cmix_s,
184 cq_s,
185 perf_rates,
186 deferred_logger);
187 }
188
189
190
191 template<typename TypeTag>
192 template<class Value>
193 void
195 computePerfRate(const std::vector<Value>& mob,
196 const Value& pressure,
197 const Value& bhp,
198 const Value& rs,
199 const Value& rv,
200 const Value& rvw,
201 const Value& rsw,
202 std::vector<Value>& b_perfcells_dense,
203 const std::vector<Scalar>& Tw,
204 const int perf,
205 const bool allow_cf,
206 const Value& skin_pressure,
207 const std::vector<Value>& cmix_s,
208 std::vector<Value>& cq_s,
209 PerforationRates<Scalar>& perf_rates,
210 DeferredLogger& deferred_logger) const
211 {
212 // Pressure drawdown (also used to determine direction of flow)
213 const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
214 Value drawdown = pressure - well_pressure;
215 if (this->isInjector()) {
216 drawdown += skin_pressure;
217 }
218
219 RatioCalculator<Value> ratioCalc{
220 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
221 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx)
222 : -1,
223 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
224 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx)
225 : -1,
226 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
227 ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx)
228 : -1,
229 this->name()
230 };
231
232 // producing perforations
233 if (drawdown > 0) {
234 // Do nothing if crossflow is not allowed
235 if (!allow_cf && this->isInjector()) {
236 return;
237 }
238
239 // compute component volumetric rates at standard conditions
240 for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
241 const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
242 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
243 }
244
245 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
246 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
247 {
248 ratioCalc.gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw,
249 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
250 this->isProducer());
251 } else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
252 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
253 {
254 ratioCalc.gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw, this->isProducer());
255 }
256 } else {
257 // Do nothing if crossflow is not allowed
258 if (!allow_cf && this->isProducer()) {
259 return;
260 }
261
262 // Using total mobilities
263 Value total_mob_dense = mob[0];
264 for (int componentIdx = 1; componentIdx < this->numConservationQuantities(); ++componentIdx) {
265 total_mob_dense += mob[componentIdx];
266 }
267
268 // compute volume ratio between connection at standard conditions
269 Value volumeRatio = bhp * 0.0; // initialize it with the correct type
270
271 if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
272 ratioCalc.disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
273 cmix_s, b_perfcells_dense, deferred_logger);
274 // DISGASW only supported for gas-water CO2STORE/H2STORE case
275 // and the simulator will throw long before it reach to this point in the code
276 // For blackoil support of DISGASW we need to add the oil component here
277 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
278 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
279 assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
280 } else {
281
282 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
283 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
284 volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
285 }
286
287 if constexpr (Indices::enableSolvent) {
288 volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
289 }
290
291 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
292 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
293 {
294 ratioCalc.gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
295 cmix_s, b_perfcells_dense,
296 deferred_logger);
297 } else {
298 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
299 const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
300 volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
301 }
302 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
303 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
304 volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
305 }
306 }
307 }
308
309 // injecting connections total volumerates at standard conditions
310 for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
311 const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
312 Value cqt_is = cqt_i / volumeRatio;
313 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
314 }
315
316 // calculating the perforation solution gas rate and solution oil rates
317 if (this->isProducer()) {
318 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
319 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
320 {
321 ratioCalc.gasOilPerfRateInj(cq_s, perf_rates,
322 rv, rs, pressure, rvw,
323 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
324 deferred_logger);
325 }
326 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
327 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
328 {
329 //no oil
330 ratioCalc.gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
331 pressure, deferred_logger);
332 }
333 }
334 }
335 }
336
337
338 template<typename TypeTag>
339 void
342 const double dt,
343 const Well::InjectionControls& inj_controls,
344 const Well::ProductionControls& prod_controls,
345 WellStateType& well_state,
346 const GroupState<Scalar>& group_state,
347 DeferredLogger& deferred_logger)
348 {
349 // TODO: only_wells should be put back to save some computation
350 // for example, the matrices B C does not need to update if only_wells
351 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
352
353 // clear all entries
354 this->linSys_.clear();
355
356 assembleWellEqWithoutIterationImpl(simulator, dt, inj_controls,
357 prod_controls, well_state,
358 group_state, deferred_logger);
359 }
360
361
362
363
364 template<typename TypeTag>
365 void
368 const double dt,
369 const Well::InjectionControls& inj_controls,
370 const Well::ProductionControls& prod_controls,
371 WellStateType& well_state,
372 const GroupState<Scalar>& group_state,
373 DeferredLogger& deferred_logger)
374 {
375 // try to regularize equation if the well does not converge
376 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
377 const Scalar volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
378
379 auto& ws = well_state.well(this->index_of_well_);
380 ws.phase_mixing_rates.fill(0.0);
381
382
383 const int np = this->number_of_phases_;
384
385 std::vector<RateVector> connectionRates = this->connectionRates_; // Copy to get right size.
386
387 auto& perf_data = ws.perf_data;
388 auto& perf_rates = perf_data.phase_rates;
389 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
390 // Calculate perforation quantities.
391 std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
392 EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
393 EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
394 calculateSinglePerf(simulator, perf, well_state, connectionRates,
395 cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
396
397 // Equation assembly for this perforation.
398 if constexpr (has_polymer && Base::has_polymermw) {
399 if (this->isInjector()) {
400 handleInjectivityEquations(simulator, well_state, perf,
401 water_flux_s, deferred_logger);
402 }
403 }
404 for (int componentIdx = 0; componentIdx < this->num_conservation_quantities_; ++componentIdx) {
405 // the cq_s entering mass balance equations need to consider the efficiency factors.
406 const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
407
408 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
409
411 assemblePerforationEq(cq_s_effective,
412 componentIdx,
413 perf,
414 this->primary_variables_.numWellEq(),
415 this->linSys_);
416
417 // Store the perforation phase flux for later usage.
418 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
419 auto& perf_rate_solvent = perf_data.solvent_rates;
420 perf_rate_solvent[perf] = cq_s[componentIdx].value();
421 } else {
422 perf_rates[perf*np + FluidSystem::activeCompToActivePhaseIdx(componentIdx)] = cq_s[componentIdx].value();
423 }
424 }
425
426 if constexpr (has_zFraction) {
428 assembleZFracEq(cq_s_zfrac_effective,
429 perf,
430 this->primary_variables_.numWellEq(),
431 this->linSys_);
432 }
433 }
434 // Update the connection
435 this->connectionRates_ = connectionRates;
436
437 // Accumulate dissolved gas and vaporized oil flow rates across all
438 // ranks sharing this well (this->index_of_well_).
439 {
440 const auto& comm = this->parallel_well_info_.communication();
441 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
442 }
443
444 // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
445 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
446
447 // add vol * dF/dt + Q to the well equations;
448 for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
449 // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
450 // since all the rates are under surface condition
451 EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
452 if (FluidSystem::numActivePhases() > 1) {
453 assert(dt > 0);
454 resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
455 this->F0_[componentIdx]) * volume / dt;
456 }
457 resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
459 assembleSourceEq(resWell_loc,
460 componentIdx,
461 this->primary_variables_.numWellEq(),
462 this->linSys_);
463 }
464
465 const auto& summaryState = simulator.vanguard().summaryState();
466 const Schedule& schedule = simulator.vanguard().schedule();
467 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
469 assembleControlEq(well_state, group_state,
470 schedule, summaryState,
471 inj_controls, prod_controls,
472 this->primary_variables_,
473 this->connections_.rho(),
474 this->linSys_,
475 stopped_or_zero_target,
476 deferred_logger);
477
478
479 // do the local inversion of D.
480 try {
481 this->linSys_.invert();
482 } catch( ... ) {
483 OPM_DEFLOG_PROBLEM(NumericalProblem, "Error when inverting local well equations for well " + name(), deferred_logger);
484 }
485 }
486
487
488
489
490 template<typename TypeTag>
491 void
493 calculateSinglePerf(const Simulator& simulator,
494 const int perf,
495 WellStateType& well_state,
496 std::vector<RateVector>& connectionRates,
497 std::vector<EvalWell>& cq_s,
498 EvalWell& water_flux_s,
499 EvalWell& cq_s_zfrac_effective,
500 DeferredLogger& deferred_logger) const
501 {
502 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
503 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
504 const int cell_idx = this->well_cells_[perf];
505 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
506 std::vector<EvalWell> mob(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
507 getMobility(simulator, perf, mob, deferred_logger);
508
509 PerforationRates<Scalar> perf_rates;
510 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
511 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
512 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
513 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
514 cq_s, perf_rates, deferred_logger);
515
516 auto& ws = well_state.well(this->index_of_well_);
517 auto& perf_data = ws.perf_data;
518 if constexpr (has_polymer && Base::has_polymermw) {
519 if (this->isInjector()) {
520 // Store the original water flux computed from the reservoir quantities.
521 // It will be required to assemble the injectivity equations.
522 const unsigned water_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
523 water_flux_s = cq_s[water_comp_idx];
524 // Modify the water flux for the rest of this function to depend directly on the
525 // local water velocity primary variable.
526 handleInjectivityRate(simulator, perf, cq_s);
527 }
528 }
529
530 // updating the solution gas rate and solution oil rate
531 if (this->isProducer()) {
532 ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
533 ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
534 ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
535 ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
536 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perf_rates.dis_gas;
537 perf_data.phase_mixing_rates[perf][ws.dissolved_gas_in_water] = perf_rates.dis_gas_in_water;
538 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perf_rates.vap_oil;
539 perf_data.phase_mixing_rates[perf][ws.vaporized_water] = perf_rates.vap_wat;
540 }
541
542 if constexpr (has_energy) {
543 connectionRates[perf][Indices::contiEnergyEqIdx] =
544 connectionRateEnergy(cq_s, intQuants, deferred_logger);
545 }
546
547 if constexpr (has_polymer) {
548 std::variant<Scalar,EvalWell> polymerConcentration;
549 if (this->isInjector()) {
550 polymerConcentration = this->wpolymer();
551 } else {
552 polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
553 intQuants.polymerViscosityCorrection());
554 }
555
556 [[maybe_unused]] EvalWell cq_s_poly;
557 std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
558 cq_s_poly) =
559 this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
560 cq_s, polymerConcentration);
561
562 if constexpr (Base::has_polymermw) {
563 updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
564 perf, connectionRates, deferred_logger);
565 }
566 }
567
568 if constexpr (has_foam) {
569 std::variant<Scalar,EvalWell> foamConcentration;
570 if (this->isInjector()) {
571 foamConcentration = this->wfoam();
572 } else {
573 foamConcentration = this->extendEval(intQuants.foamConcentration());
574 }
575 connectionRates[perf][Indices::contiFoamEqIdx] =
576 this->connections_.connectionRateFoam(cq_s, foamConcentration,
577 FoamModule::transportPhase(),
578 deferred_logger);
579 }
580
581 if constexpr (has_zFraction) {
582 std::variant<Scalar,std::array<EvalWell,2>> solventConcentration;
583 if (this->isInjector()) {
584 solventConcentration = this->wsolvent();
585 } else {
586 solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
587 this->extendEval(intQuants.yVolume())};
588 }
589 std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
590 cq_s_zfrac_effective) =
591 this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
592 perf_rates.dis_gas, cq_s,
593 solventConcentration);
594 }
595
596 if constexpr (has_brine) {
597 std::variant<Scalar,EvalWell> saltConcentration;
598 if (this->isInjector()) {
599 saltConcentration = this->wsalt();
600 } else {
601 saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
602 }
603
604 connectionRates[perf][Indices::contiBrineEqIdx] =
605 this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
606 perf_rates.vap_wat, cq_s,
607 saltConcentration);
608 }
609
610 if constexpr (has_micp) {
611 std::variant<Scalar,EvalWell> microbialConcentration;
612 std::variant<Scalar,EvalWell> oxygenConcentration;
613 std::variant<Scalar,EvalWell> ureaConcentration;
614 if (this->isInjector()) {
615 microbialConcentration = this->wmicrobes();
616 oxygenConcentration = this->woxygen();
617 ureaConcentration = this->wurea();
618 } else {
619 microbialConcentration = this->extendEval(intQuants.microbialConcentration());
620 oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
621 ureaConcentration = this->extendEval(intQuants.ureaConcentration());
622 }
623 std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
624 connectionRates[perf][Indices::contiOxygenEqIdx],
625 connectionRates[perf][Indices::contiUreaEqIdx]) =
626 this->connections_.connectionRatesMICP(perf_data.microbial_rates[perf],
627 perf_data.oxygen_rates[perf],
628 perf_data.urea_rates[perf],
629 cq_s,
630 microbialConcentration,
631 oxygenConcentration,
632 ureaConcentration);
633 }
634
635 // Store the perforation pressure for later usage.
636 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
637
638 // Store the perforation gass mass rate.
639 if (FluidSystem::phaseUsage().hasCO2orH2Store()) {
640 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
641 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
642 perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
643 }
644
645 // Store the perforation water mass rate.
646 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
647 const unsigned wat_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
648 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::waterPhaseIdx, Base::pvtRegionIdx() );
649 perf_data.wat_mass_rates[perf] = cq_s[wat_comp_idx].value() * rho;
650 }
651 }
652
653
654
655 template<typename TypeTag>
656 template<class Value>
657 void
659 getMobility(const Simulator& simulator,
660 const int perf,
661 std::vector<Value>& mob,
662 DeferredLogger& deferred_logger) const
663 {
664 auto obtain = [this](const Eval& value)
665 {
666 if constexpr (std::is_same_v<Value, Scalar>) {
667 static_cast<void>(this); // suppress clang warning
668 return getValue(value);
669 } else {
670 return this->extendEval(value);
671 }
672 };
673 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
674 obtain, deferred_logger);
675
676 // modify the water mobility if polymer is present
677 if constexpr (has_polymer) {
678 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
679 OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
680 }
681
682 // for the cases related to polymer molecular weight, we assume fully mixing
683 // as a result, the polymer and water share the same viscosity
684 if constexpr (!Base::has_polymermw) {
685 if constexpr (std::is_same_v<Value, Scalar>) {
686 std::vector<EvalWell> mob_eval(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
687 for (std::size_t i = 0; i < mob.size(); ++i) {
688 mob_eval[i].setValue(mob[i]);
689 }
690 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
691 for (std::size_t i = 0; i < mob.size(); ++i) {
692 mob[i] = getValue(mob_eval[i]);
693 }
694 } else {
695 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
696 }
697 }
698 }
699
700 // if the injecting well has WINJMULT setup, we update the mobility accordingly
701 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
702 const Scalar bhp = this->primary_variables_.value(Bhp);
703 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
704 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
705 for (std::size_t i = 0; i < mob.size(); ++i) {
706 mob[i] *= multiplier;
707 }
708 }
709 }
710
711
712 template<typename TypeTag>
713 void
715 updateWellState(const Simulator& simulator,
716 const BVectorWell& dwells,
717 WellStateType& well_state,
718 DeferredLogger& deferred_logger)
719 {
720 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
721
722 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
723 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
724
725 const auto& summary_state = simulator.vanguard().summaryState();
726 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
727 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
728 }
729
730
731
732
733
734 template<typename TypeTag>
735 void
738 const bool stop_or_zero_rate_target,
739 DeferredLogger& deferred_logger)
740 {
741 const Scalar dFLimit = this->param_.dwell_fraction_max_;
742 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
743 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
744
745 // for the water velocity and skin pressure
746 if constexpr (Base::has_polymermw) {
747 this->primary_variables_.updateNewtonPolyMW(dwells);
748 }
749
750 this->primary_variables_.checkFinite(deferred_logger);
751 }
752
753
754
755
756
757 template<typename TypeTag>
758 void
761 const SummaryState& summary_state,
762 DeferredLogger& deferred_logger) const
763 {
764 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
765
766 // other primary variables related to polymer injectivity study
767 if constexpr (Base::has_polymermw) {
768 this->primary_variables_.copyToWellStatePolyMW(well_state);
769 }
770 }
771
772
773
774
775
776 template<typename TypeTag>
777 void
779 updateIPR(const Simulator& simulator, DeferredLogger& deferred_logger) const
780 {
781 // TODO: not handling solvent related here for now
782
783 // initialize all the values to be zero to begin with
784 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
785 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
786
787 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
788 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
789 getMobility(simulator, perf, mob, deferred_logger);
790
791 const int cell_idx = this->well_cells_[perf];
792 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
793 const auto& fs = int_quantities.fluidState();
794 // the pressure of the reservoir grid block the well connection is in
795 Scalar p_r = this->getPerfCellPressure(fs).value();
796
797 // calculating the b for the connection
798 std::vector<Scalar> b_perf(this->num_conservation_quantities_);
799 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
800 if (!FluidSystem::phaseIsActive(phase)) {
801 continue;
802 }
803 const unsigned comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phase));
804 b_perf[comp_idx] = fs.invB(phase).value();
805 }
806 if constexpr (has_solvent) {
807 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
808 }
809
810 // the pressure difference between the connection and BHP
811 const Scalar h_perf = this->connections_.pressure_diff(perf);
812 const Scalar pressure_diff = p_r - h_perf;
813
814 // Let us add a check, since the pressure is calculated based on zero value BHP
815 // it should not be negative anyway. If it is negative, we might need to re-formulate
816 // to taking into consideration the crossflow here.
817 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
818 deferred_logger.debug("CROSSFLOW_IPR",
819 "cross flow found when updateIPR for well " + name()
820 + " . The connection is ignored in IPR calculations");
821 // we ignore these connections for now
822 continue;
823 }
824
825 // the well index associated with the connection
826 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
827 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
828 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
829 int_quantities,
830 trans_mult,
831 wellstate_nupcol);
832 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
833 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
834 for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
835 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
836 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
837 ipr_b_perf[comp_idx] += tw_mob;
838 }
839
840 // we need to handle the rs and rv when both oil and gas are present
841 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
842 const unsigned oil_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
843 const unsigned gas_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
844 const Scalar rs = (fs.Rs()).value();
845 const Scalar rv = (fs.Rv()).value();
846
847 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
848 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
849
850 ipr_a_perf[gas_comp_idx] += dis_gas_a;
851 ipr_a_perf[oil_comp_idx] += vap_oil_a;
852
853 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
854 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
855
856 ipr_b_perf[gas_comp_idx] += dis_gas_b;
857 ipr_b_perf[oil_comp_idx] += vap_oil_b;
858 }
859
860 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
861 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
862 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
863 }
864 }
865 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
866 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
867 }
868
869 template<typename TypeTag>
870 void
872 updateIPRImplicit(const Simulator& simulator,
873 WellStateType& well_state,
874 DeferredLogger& deferred_logger)
875 {
876 // Compute IPR based on *converged* well-equation:
877 // For a component rate r the derivative dr/dbhp is obtained by
878 // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target)
879 // where Eq(x)=0 is the well equation setup with bhp control and primary variables x
880
881 // We shouldn't have zero rates at this stage, but check
882 bool zero_rates;
883 auto rates = well_state.well(this->index_of_well_).surface_rates;
884 zero_rates = true;
885 for (std::size_t p = 0; p < rates.size(); ++p) {
886 zero_rates &= rates[p] == 0.0;
887 }
888 auto& ws = well_state.well(this->index_of_well_);
889 if (zero_rates) {
890 const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
891 deferred_logger.debug(msg);
892 /*
893 // could revert to standard approach here:
894 updateIPR(simulator, deferred_logger);
895 for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
896 const int idx = this->activeCompToActivePhaseIdx(comp_idx);
897 ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
898 ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
899 }
900 return;
901 */
902 }
903 const auto& group_state = simulator.problem().wellModel().groupState();
904
905 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
906 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
907
908 auto inj_controls = Well::InjectionControls(0);
909 auto prod_controls = Well::ProductionControls(0);
910 prod_controls.addControl(Well::ProducerCMode::BHP);
911 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
912
913 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
914 const auto cmode = ws.production_cmode;
915 ws.production_cmode = Well::ProducerCMode::BHP;
916 const double dt = simulator.timeStepSize();
917 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
918
919 const size_t nEq = this->primary_variables_.numWellEq();
920 BVectorWell rhs(1);
921 rhs[0].resize(nEq);
922 // rhs = 0 except -1 for control eq
923 for (size_t i=0; i < nEq; ++i){
924 rhs[0][i] = 0.0;
925 }
926 rhs[0][Bhp] = -1.0;
927
928 BVectorWell x_well(1);
929 x_well[0].resize(nEq);
930 this->linSys_.solve(rhs, x_well);
931
932 for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
933 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
934 const int idx = FluidSystem::activeCompToActivePhaseIdx(comp_idx);
935 for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
936 // well primary variable derivatives in EvalWell start at position Indices::numEq
937 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
938 }
939 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
940 }
941 // reset cmode
942 ws.production_cmode = cmode;
943 }
944
945 template<typename TypeTag>
946 void
949 const Simulator& simulator,
950 DeferredLogger& deferred_logger)
951 {
952 const auto& summaryState = simulator.vanguard().summaryState();
953 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
954 // Crude but works: default is one atmosphere.
955 // TODO: a better way to detect whether the BHP is defaulted or not
956 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
957 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
958 // if the BHP limit is not defaulted or the well does not have a THP limit
959 // we need to check the BHP limit
960 Scalar total_ipr_mass_rate = 0.0;
961 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
962 {
963 if (!FluidSystem::phaseIsActive(phaseIdx)) {
964 continue;
965 }
966
967 const unsigned compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
968 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
969
970 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
971 total_ipr_mass_rate += ipr_rate * rho;
972 }
973 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
974 this->operability_status_.operable_under_only_bhp_limit = false;
975 }
976
977 // checking whether running under BHP limit will violate THP limit
978 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
979 // option 1: calculate well rates based on the BHP limit.
980 // option 2: stick with the above IPR curve
981 // we use IPR here
982 std::vector<Scalar> well_rates_bhp_limit;
983 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
984
985 this->adaptRatesForVFP(well_rates_bhp_limit);
986 const Scalar thp_limit = this->getTHPConstraint(summaryState);
987 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
988 bhp_limit,
989 this->connections_.rho(),
990 this->getALQ(well_state),
991 thp_limit,
992 deferred_logger);
993 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
994 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
995 }
996 }
997 } else {
998 // defaulted BHP and there is a THP constraint
999 // default BHP limit is about 1 atm.
1000 // when applied the hydrostatic pressure correction dp,
1001 // most likely we get a negative value (bhp + dp)to search in the VFP table,
1002 // which is not desirable.
1003 // we assume we can operate under defaulted BHP limit and will violate the THP limit
1004 // when operating under defaulted BHP limit.
1005 this->operability_status_.operable_under_only_bhp_limit = true;
1006 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1007 }
1008 }
1009
1010
1011
1012
1013
1014 template<typename TypeTag>
1015 void
1018 const WellStateType& well_state,
1019 DeferredLogger& deferred_logger)
1020 {
1021 const auto& summaryState = simulator.vanguard().summaryState();
1022 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1023 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1024
1025 if (obtain_bhp) {
1026 this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1027
1028 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1029 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1030 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1031
1032 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1033 if (this->isProducer() && *obtain_bhp < thp_limit) {
1034 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1035 + " bars is SMALLER than thp limit "
1036 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1037 + " bars as a producer for well " + name();
1038 deferred_logger.debug(msg);
1039 }
1040 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1041 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1042 + " bars is LARGER than thp limit "
1043 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1044 + " bars as a injector for well " + name();
1045 deferred_logger.debug(msg);
1046 }
1047 } else {
1048 this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1049 this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1050 if (!this->wellIsStopped()) {
1051 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1052 deferred_logger.debug(" could not find bhp value at thp limit "
1053 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1054 + " bar for well " + name() + ", the well might need to be closed ");
1055 }
1056 }
1057 }
1058
1059
1060
1061
1062
1063 template<typename TypeTag>
1064 bool
1066 allDrawDownWrongDirection(const Simulator& simulator) const
1067 {
1068 bool all_drawdown_wrong_direction = true;
1069
1070 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1071 const int cell_idx = this->well_cells_[perf];
1072 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1073 const auto& fs = intQuants.fluidState();
1074
1075 const Scalar pressure = this->getPerfCellPressure(fs).value();
1076 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1077
1078 // Pressure drawdown (also used to determine direction of flow)
1079 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1080 const Scalar drawdown = pressure - well_pressure;
1081
1082 // for now, if there is one perforation can produce/inject in the correct
1083 // direction, we consider this well can still produce/inject.
1084 // TODO: it can be more complicated than this to cause wrong-signed rates
1085 if ( (drawdown < 0. && this->isInjector()) ||
1086 (drawdown > 0. && this->isProducer()) ) {
1087 all_drawdown_wrong_direction = false;
1088 break;
1089 }
1090 }
1091
1092 const auto& comm = this->parallel_well_info_.communication();
1093 if (comm.size() > 1)
1094 {
1095 all_drawdown_wrong_direction =
1096 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1097 }
1098
1099 return all_drawdown_wrong_direction;
1100 }
1101
1102
1103
1104
1105 template<typename TypeTag>
1106 bool
1108 openCrossFlowAvoidSingularity(const Simulator& simulator) const
1109 {
1110 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1111 }
1112
1113
1114
1115
1116 template<typename TypeTag>
1120 const WellStateType& well_state) const
1121 {
1122 auto prop_func = typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1123 // getTemperature
1124 [&model = simulator.model()](int cell_idx, int phase_idx)
1125 {
1126 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1127 .fluidState().temperature(phase_idx).value();
1128 },
1129
1130 // getSaltConcentration
1131 [&model = simulator.model()](int cell_idx)
1132 {
1133 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1134 .fluidState().saltConcentration().value();
1135 },
1136
1137 // getPvtRegionIdx
1138 [&model = simulator.model()](int cell_idx)
1139 {
1140 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1141 .fluidState().pvtRegionIndex();
1142 }
1143 };
1144
1145 if constexpr (Indices::enableSolvent) {
1146 prop_func.solventInverseFormationVolumeFactor =
1147 [&model = simulator.model()](int cell_idx)
1148 {
1149 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1150 .solventInverseFormationVolumeFactor().value();
1151 };
1152
1153 prop_func.solventRefDensity = [&model = simulator.model()](int cell_idx)
1154 {
1155 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1156 .solventRefDensity();
1157 };
1158 }
1159
1160 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1161 }
1162
1163
1164
1165
1166
1167 template<typename TypeTag>
1170 getWellConvergence(const Simulator& simulator,
1171 const WellStateType& well_state,
1172 const std::vector<Scalar>& B_avg,
1173 DeferredLogger& deferred_logger,
1174 const bool relax_tolerance) const
1175 {
1176 // the following implementation assume that the polymer is always after the w-o-g phases
1177 // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
1178 assert((int(B_avg.size()) == this->num_conservation_quantities_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1179
1180 Scalar tol_wells = this->param_.tolerance_wells_;
1181 // use stricter tolerance for stopped wells and wells under zero rate target control.
1182 constexpr Scalar stopped_factor = 1.e-4;
1183 // use stricter tolerance for dynamic thp to ameliorate network convergence
1184 constexpr Scalar dynamic_thp_factor = 1.e-1;
1185 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1186 tol_wells = tol_wells*stopped_factor;
1187 } else if (this->getDynamicThpLimit()) {
1188 tol_wells = tol_wells*dynamic_thp_factor;
1189 }
1190
1191 std::vector<Scalar> res;
1192 ConvergenceReport report = this->StdWellEval::getWellConvergence(well_state,
1193 B_avg,
1194 this->param_.max_residual_allowed_,
1195 tol_wells,
1196 this->param_.relaxed_tolerance_flow_well_,
1197 relax_tolerance,
1198 this->wellIsStopped(),
1199 res,
1200 deferred_logger);
1201
1202 checkConvergenceExtraEqs(res, report);
1203
1204 return report;
1205 }
1206
1207
1208
1209
1210
1211 template<typename TypeTag>
1212 void
1214 updateProductivityIndex(const Simulator& simulator,
1215 const WellProdIndexCalculator<Scalar>& wellPICalc,
1216 WellStateType& well_state,
1217 DeferredLogger& deferred_logger) const
1218 {
1219 auto fluidState = [&simulator, this](const int perf)
1220 {
1221 const auto cell_idx = this->well_cells_[perf];
1222 return simulator.model()
1223 .intensiveQuantities(cell_idx, /*timeIdx=*/ 0).fluidState();
1224 };
1225
1226 const int np = this->number_of_phases_;
1227 auto setToZero = [np](Scalar* x) -> void
1228 {
1229 std::fill_n(x, np, 0.0);
1230 };
1231
1232 auto addVector = [np](const Scalar* src, Scalar* dest) -> void
1233 {
1234 std::transform(src, src + np, dest, dest, std::plus<>{});
1235 };
1236
1237 auto& ws = well_state.well(this->index_of_well_);
1238 auto& perf_data = ws.perf_data;
1239 auto* wellPI = ws.productivity_index.data();
1240 auto* connPI = perf_data.prod_index.data();
1241
1242 setToZero(wellPI);
1243
1244 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1245 auto subsetPerfID = 0;
1246
1247 for (const auto& perf : *this->perf_data_) {
1248 auto allPerfID = perf.ecl_index;
1249
1250 auto connPICalc = [&wellPICalc, allPerfID](const Scalar mobility) -> Scalar
1251 {
1252 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
1253 };
1254
1255 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
1256 getMobility(simulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
1257
1258 const auto& fs = fluidState(subsetPerfID);
1259 setToZero(connPI);
1260
1261 if (this->isInjector()) {
1262 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1263 mob, connPI, deferred_logger);
1264 }
1265 else { // Production or zero flow rate
1266 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1267 }
1268
1269 addVector(connPI, wellPI);
1270
1271 ++subsetPerfID;
1272 connPI += np;
1273 }
1274
1275 // Sum with communication in case of distributed well.
1276 const auto& comm = this->parallel_well_info_.communication();
1277 if (comm.size() > 1) {
1278 comm.sum(wellPI, np);
1279 }
1280
1281 assert ((static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1282 "Internal logic error in processing connections for PI/II");
1283 }
1284
1285
1286
1287 template<typename TypeTag>
1290 const WellStateType& well_state,
1291 const WellConnectionProps& props,
1292 DeferredLogger& deferred_logger)
1293 {
1294 // Cell level dynamic property call-back functions as fall-back
1295 // option for calculating connection level mixture densities in
1296 // stopped or zero-rate producer wells.
1297 const auto prop_func = typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1298 // This becomes slightly more palatable with C++20's designated
1299 // initialisers.
1300
1301 // mobility: Phase mobilities in specified cell.
1302 [&model = simulator.model()](const int cell,
1303 const std::vector<int>& phases,
1304 std::vector<Scalar>& mob)
1305 {
1306 const auto& iq = model.intensiveQuantities(cell, /* time_idx = */ 0);
1307
1308 std::transform(phases.begin(), phases.end(), mob.begin(),
1309 [&iq](const int phase) { return iq.mobility(phase).value(); });
1310 },
1311
1312 // densityInCell: Reservoir condition phase densities in
1313 // specified cell.
1314 [&model = simulator.model()](const int cell,
1315 const std::vector<int>& phases,
1316 std::vector<Scalar>& rho)
1317 {
1318 const auto& fs = model.intensiveQuantities(cell, /* time_idx = */ 0).fluidState();
1319
1320 std::transform(phases.begin(), phases.end(), rho.begin(),
1321 [&fs](const int phase) { return fs.density(phase).value(); });
1322 }
1323 };
1324
1325 const auto stopped_or_zero_rate_target = this->
1326 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1327
1328 this->connections_
1329 .computeProperties(stopped_or_zero_rate_target, well_state,
1330 prop_func, props, deferred_logger);
1331 }
1332
1333
1334
1335
1336
1337 template<typename TypeTag>
1338 void
1341 const WellStateType& well_state,
1342 DeferredLogger& deferred_logger)
1343 {
1344 const auto props = computePropertiesForWellConnectionPressures
1345 (simulator, well_state);
1346
1347 computeWellConnectionDensitesPressures(simulator, well_state,
1348 props, deferred_logger);
1349 }
1350
1351
1352
1353
1354
1355 template<typename TypeTag>
1356 void
1358 solveEqAndUpdateWellState(const Simulator& simulator,
1359 WellStateType& well_state,
1360 DeferredLogger& deferred_logger)
1361 {
1362 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1363
1364 // We assemble the well equations, then we check the convergence,
1365 // which is why we do not put the assembleWellEq here.
1366 BVectorWell dx_well(1);
1367 dx_well[0].resize(this->primary_variables_.numWellEq());
1368 this->linSys_.solve( dx_well);
1369
1370 updateWellState(simulator, dx_well, well_state, deferred_logger);
1371 }
1372
1373
1374
1375
1376
1377 template<typename TypeTag>
1378 void
1381 const WellStateType& well_state,
1382 DeferredLogger& deferred_logger)
1383 {
1384 updatePrimaryVariables(simulator, well_state, deferred_logger);
1385 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1386 this->computeAccumWell();
1387 }
1388
1389
1390
1391 template<typename TypeTag>
1392 void
1394 apply(const BVector& x, BVector& Ax) const
1395 {
1396 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1397
1398 if (this->param_.matrix_add_well_contributions_)
1399 {
1400 // Contributions are already in the matrix itself
1401 return;
1402 }
1403
1404 this->linSys_.apply(x, Ax);
1405 }
1406
1407
1408
1409
1410 template<typename TypeTag>
1411 void
1413 apply(BVector& r) const
1414 {
1415 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1416
1417 this->linSys_.apply(r);
1418 }
1419
1420
1421
1422
1423 template<typename TypeTag>
1424 void
1427 const BVector& x,
1428 WellStateType& well_state,
1429 DeferredLogger& deferred_logger)
1430 {
1431 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1432
1433 BVectorWell xw(1);
1434 xw[0].resize(this->primary_variables_.numWellEq());
1435
1436 this->linSys_.recoverSolutionWell(x, xw);
1437 updateWellState(simulator, xw, well_state, deferred_logger);
1438 }
1439
1440
1441
1442
1443 template<typename TypeTag>
1444 void
1446 computeWellRatesWithBhp(const Simulator& simulator,
1447 const Scalar& bhp,
1448 std::vector<Scalar>& well_flux,
1449 DeferredLogger& deferred_logger) const
1450 {
1451 OPM_TIMEFUNCTION();
1452 const int np = this->number_of_phases_;
1453 well_flux.resize(np, 0.0);
1454
1455 const bool allow_cf = this->getAllowCrossFlow();
1456
1457 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1458 const int cell_idx = this->well_cells_[perf];
1459 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1460 // flux for each perforation
1461 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
1462 getMobility(simulator, perf, mob, deferred_logger);
1463 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1464 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1465 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1466
1467 std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
1468 PerforationRates<Scalar> perf_rates;
1469 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1470 cq_s, perf_rates, deferred_logger);
1471
1472 for(int p = 0; p < np; ++p) {
1473 well_flux[FluidSystem::activeCompToActivePhaseIdx(p)] += cq_s[p];
1474 }
1475
1476 // the solvent contribution is added to the gas potentials
1477 if constexpr (has_solvent) {
1478 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
1479 // TODO: should we use compIdx here?
1480 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1481 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1482 }
1483 }
1484 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1485 }
1486
1487
1488
1489 template<typename TypeTag>
1490 void
1493 const Scalar& bhp,
1494 std::vector<Scalar>& well_flux,
1495 DeferredLogger& deferred_logger) const
1496 {
1497 // creating a copy of the well itself, to avoid messing up the explicit information
1498 // during this copy, the only information not copied properly is the well controls
1499 StandardWell<TypeTag> well_copy(*this);
1500 well_copy.resetDampening();
1501
1502 // iterate to get a more accurate well density
1503 // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
1504 // to replace the original one
1505 WellStateType well_state_copy = simulator.problem().wellModel().wellState();
1506 const auto& group_state = simulator.problem().wellModel().groupState();
1507
1508 // Get the current controls.
1509 const auto& summary_state = simulator.vanguard().summaryState();
1510 auto inj_controls = well_copy.well_ecl_.isInjector()
1511 ? well_copy.well_ecl_.injectionControls(summary_state)
1512 : Well::InjectionControls(0);
1513 auto prod_controls = well_copy.well_ecl_.isProducer()
1514 ? well_copy.well_ecl_.productionControls(summary_state) :
1515 Well::ProductionControls(0);
1516
1517 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
1518 auto& ws = well_state_copy.well(this->index_of_well_);
1519 if (well_copy.well_ecl_.isInjector()) {
1520 inj_controls.bhp_limit = bhp;
1521 ws.injection_cmode = Well::InjectorCMode::BHP;
1522 } else {
1523 prod_controls.bhp_limit = bhp;
1524 ws.production_cmode = Well::ProducerCMode::BHP;
1525 }
1526 ws.bhp = bhp;
1527
1528 // initialized the well rates with the potentials i.e. the well rates based on bhp
1529 const int np = this->number_of_phases_;
1530 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1531 for (int phase = 0; phase < np; ++phase){
1532 well_state_copy.wellRates(this->index_of_well_)[phase]
1533 = sign * ws.well_potentials[phase];
1534 }
1535 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1536 well_copy.computeAccumWell();
1537
1538 const double dt = simulator.timeStepSize();
1539 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1540 if (!converged) {
1541 const std::string msg = " well " + name() + " did not get converged during well potential calculations "
1542 " potentials are computed based on unconverged solution";
1543 deferred_logger.debug(msg);
1544 }
1545 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1546 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1547 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1548 }
1549
1550
1551
1552
1553 template<typename TypeTag>
1554 std::vector<typename StandardWell<TypeTag>::Scalar>
1557 DeferredLogger& deferred_logger,
1558 const WellStateType& well_state) const
1559 {
1560 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1561 const auto& summary_state = simulator.vanguard().summaryState();
1562
1563 const auto& well = this->well_ecl_;
1564 if (well.isInjector()){
1565 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1566 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1567 if (bhp_at_thp_limit) {
1568 const Scalar bhp = std::min(*bhp_at_thp_limit,
1569 static_cast<Scalar>(controls.bhp_limit));
1570 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1571 } else {
1572 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1573 "Failed in getting converged thp based potential calculation for well "
1574 + name() + ". Instead the bhp based value is used");
1575 const Scalar bhp = controls.bhp_limit;
1576 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1577 }
1578 } else {
1579 computeWellRatesWithThpAlqProd(
1580 simulator, summary_state,
1581 deferred_logger, potentials, this->getALQ(well_state)
1582 );
1583 }
1584
1585 return potentials;
1586 }
1587
1588 template<typename TypeTag>
1589 bool
1592 const WellStateType& well_state,
1593 std::vector<Scalar>& well_potentials,
1594 DeferredLogger& deferred_logger) const
1595 {
1596 // Create a copy of the well.
1597 // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials
1598 // is allready a copy, but not from other calls.
1599 StandardWell<TypeTag> well_copy(*this);
1600
1601 // store a copy of the well state, we don't want to update the real well state
1602 WellStateType well_state_copy = well_state;
1603 const auto& group_state = simulator.problem().wellModel().groupState();
1604 auto& ws = well_state_copy.well(this->index_of_well_);
1605
1606 // get current controls
1607 const auto& summary_state = simulator.vanguard().summaryState();
1608 auto inj_controls = well_copy.well_ecl_.isInjector()
1609 ? well_copy.well_ecl_.injectionControls(summary_state)
1610 : Well::InjectionControls(0);
1611 auto prod_controls = well_copy.well_ecl_.isProducer()
1612 ? well_copy.well_ecl_.productionControls(summary_state) :
1613 Well::ProductionControls(0);
1614
1615 // prepare/modify well state and control
1616 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1617
1618 // update connection pressures relative to updated bhp to get better estimate of connection dp
1619 const int num_perf = ws.perf_data.size();
1620 for (int perf = 0; perf < num_perf; ++perf) {
1621 ws.perf_data.pressure[perf] = ws.bhp + well_copy.connections_.pressure_diff(perf);
1622 }
1623 // initialize rates from previous potentials
1624 const int np = this->number_of_phases_;
1625 bool trivial = true;
1626 for (int phase = 0; phase < np; ++phase){
1627 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1628 }
1629 if (!trivial) {
1630 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1631 for (int phase = 0; phase < np; ++phase) {
1632 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1633 }
1634 }
1635
1636 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1637 const double dt = simulator.timeStepSize();
1638 // iterate to get a solution at the given bhp.
1639 bool converged = false;
1640 if (this->well_ecl_.isProducer()) {
1641 converged = well_copy.solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1642 } else {
1643 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1644 }
1645
1646 // fetch potentials (sign is updated on the outside).
1647 well_potentials.clear();
1648 well_potentials.resize(np, 0.0);
1649 for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
1650 if (has_solvent && comp_idx == Indices::contiSolventEqIdx) continue; // we do not store the solvent in the well_potentials
1651 const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
1652 well_potentials[FluidSystem::activeCompToActivePhaseIdx(comp_idx)] = rate.value();
1653 }
1654
1655 // the solvent contribution is added to the gas potentials
1656 if constexpr (has_solvent) {
1657 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
1658 // TODO: should we use compIdx here?
1659 const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
1660 const EvalWell rate = well_copy.primary_variables_.getQs(Indices::contiSolventEqIdx);
1661 well_potentials[gas_pos] += rate.value();
1662 }
1663 return converged;
1664 }
1665
1666
1667 template<typename TypeTag>
1671 const SummaryState &summary_state,
1672 DeferredLogger& deferred_logger,
1673 std::vector<Scalar>& potentials,
1674 Scalar alq) const
1675 {
1676 Scalar bhp;
1677 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1678 simulator, summary_state, alq, deferred_logger, /*iterate_if_no_solution */ true);
1679 if (bhp_at_thp_limit) {
1680 const auto& controls = this->well_ecl_.productionControls(summary_state);
1681 bhp = std::max(*bhp_at_thp_limit,
1682 static_cast<Scalar>(controls.bhp_limit));
1683 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1684 }
1685 else {
1686 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1687 "Failed in getting converged thp based potential calculation for well "
1688 + name() + ". Instead the bhp based value is used");
1689 const auto& controls = this->well_ecl_.productionControls(summary_state);
1690 bhp = controls.bhp_limit;
1691 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1692 }
1693 return bhp;
1694 }
1695
1696 template<typename TypeTag>
1697 void
1700 const SummaryState& summary_state,
1701 DeferredLogger& deferred_logger,
1702 std::vector<Scalar>& potentials,
1703 Scalar alq) const
1704 {
1705 /*double bhp =*/
1706 computeWellRatesAndBhpWithThpAlqProd(simulator,
1707 summary_state,
1708 deferred_logger,
1709 potentials,
1710 alq);
1711 }
1712
1713 template<typename TypeTag>
1714 void
1716 computeWellPotentials(const Simulator& simulator,
1717 const WellStateType& well_state,
1718 std::vector<Scalar>& well_potentials,
1719 DeferredLogger& deferred_logger) // const
1720 {
1721 const auto [compute_potential, bhp_controlled_well] =
1723
1724 if (!compute_potential) {
1725 return;
1726 }
1727
1728 bool converged_implicit = false;
1729 // for newly opened wells we dont compute the potentials implicit
1730 // group controlled wells with defaulted guiderates will have zero targets as
1731 // the potentials are used to compute the well fractions.
1732 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1733 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1734 }
1735 if (!converged_implicit) {
1736 // does the well have a THP related constraint?
1737 const auto& summaryState = simulator.vanguard().summaryState();
1738 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1739 // get the bhp value based on the bhp constraints
1740 Scalar bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1741
1742 // In some very special cases the bhp pressure target are
1743 // temporary violated. This may lead to too small or negative potentials
1744 // that could lead to premature shutting of wells.
1745 // As a remedy the bhp that gives the largest potential is used.
1746 // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp,
1747 // and the potentials will be computed using the limit as expected.
1748 const auto& ws = well_state.well(this->index_of_well_);
1749 if (this->isInjector())
1750 bhp = std::max(ws.bhp, bhp);
1751 else
1752 bhp = std::min(ws.bhp, bhp);
1753
1754 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1755 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1756 } else {
1757 // the well has a THP related constraint
1758 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1759 }
1760 }
1761
1762 this->checkNegativeWellPotentials(well_potentials,
1763 this->param_.check_well_operability_,
1764 deferred_logger);
1765 }
1766
1767
1768
1769
1770
1771
1772
1773 template<typename TypeTag>
1776 connectionDensity([[maybe_unused]] const int globalConnIdx,
1777 const int openConnIdx) const
1778 {
1779 return (openConnIdx < 0)
1780 ? 0.0
1781 : this->connections_.rho(openConnIdx);
1782 }
1783
1784
1785
1786
1787
1788 template<typename TypeTag>
1789 void
1791 updatePrimaryVariables(const Simulator& simulator,
1792 const WellStateType& well_state,
1793 DeferredLogger& deferred_logger)
1794 {
1795 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1796
1797 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1798 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1799
1800 // other primary variables related to polymer injection
1801 if constexpr (Base::has_polymermw) {
1802 this->primary_variables_.updatePolyMW(well_state);
1803 }
1804
1805 this->primary_variables_.checkFinite(deferred_logger);
1806 }
1807
1808
1809
1810
1811 template<typename TypeTag>
1814 getRefDensity() const
1815 {
1816 return this->connections_.rho();
1817 }
1818
1819
1820
1821
1822 template<typename TypeTag>
1823 void
1826 const int perf,
1827 std::vector<EvalWell>& mob,
1828 DeferredLogger& deferred_logger) const
1829 {
1830 const int cell_idx = this->well_cells_[perf];
1831 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1832 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1833
1834 // TODO: not sure should based on the well type or injecting/producing peforations
1835 // it can be different for crossflow
1836 if (this->isInjector()) {
1837 // assume fully mixing within injecting wellbore
1838 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1839 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
1840 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) );
1841 }
1842
1843 if (PolymerModule::hasPlyshlog()) {
1844 // we do not calculate the shear effects for injection wells when they do not
1845 // inject polymer.
1846 if (this->isInjector() && this->wpolymer() == 0.) {
1847 return;
1848 }
1849 // compute the well water velocity with out shear effects.
1850 // TODO: do we need to turn on crossflow here?
1851 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1852 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1853
1854 std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1855 PerforationRates<Scalar> perf_rates;
1856 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1857 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1858 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1859 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1860 perf_rates, deferred_logger);
1861 // TODO: make area a member
1862 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1863 const auto& material_law_manager = simulator.problem().materialLawManager();
1864 const auto& scaled_drainage_info =
1865 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1866 const Scalar swcr = scaled_drainage_info.Swcr;
1867 const EvalWell poro = this->extendEval(int_quant.porosity());
1868 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1869 // guard against zero porosity and no water
1870 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1871 const unsigned waterCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
1872 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1873
1874 if (PolymerModule::hasShrate()) {
1875 // the equation for the water velocity conversion for the wells and reservoir are from different version
1876 // of implementation. It can be changed to be more consistent when possible.
1877 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1878 }
1879 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1880 int_quant.pvtRegionIndex(),
1881 water_velocity);
1882 // modify the mobility with the shear factor.
1883 mob[waterCompIdx] /= shear_factor;
1884 }
1885 }
1886
1887 template<typename TypeTag>
1888 void
1890 {
1891 this->linSys_.extract(jacobian);
1892 }
1893
1894
1895 template <typename TypeTag>
1896 void
1898 const BVector& weights,
1899 const int pressureVarIndex,
1900 const bool use_well_weights,
1901 const WellStateType& well_state) const
1902 {
1903 this->linSys_.extractCPRPressureMatrix(jacobian,
1904 weights,
1905 pressureVarIndex,
1906 use_well_weights,
1907 *this,
1908 Bhp,
1909 well_state);
1910 }
1911
1912
1913
1914 template<typename TypeTag>
1917 pskinwater(const Scalar throughput,
1918 const EvalWell& water_velocity,
1919 DeferredLogger& deferred_logger) const
1920 {
1921 if constexpr (Base::has_polymermw) {
1922 const int water_table_id = this->polymerWaterTable_();
1923 if (water_table_id <= 0) {
1924 OPM_DEFLOG_THROW(std::runtime_error,
1925 fmt::format("Unused SKPRWAT table id used for well {}", name()),
1926 deferred_logger);
1927 }
1928 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1929 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1930 // the skin pressure when injecting water, which also means the polymer concentration is zero
1931 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1932 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1933 return pskin_water;
1934 } else {
1935 OPM_DEFLOG_THROW(std::runtime_error,
1936 fmt::format("Polymermw is not activated, while injecting "
1937 "skin pressure is requested for well {}", name()),
1938 deferred_logger);
1939 }
1940 }
1941
1942
1943
1944
1945
1946 template<typename TypeTag>
1949 pskin(const Scalar throughput,
1950 const EvalWell& water_velocity,
1951 const EvalWell& poly_inj_conc,
1952 DeferredLogger& deferred_logger) const
1953 {
1954 if constexpr (Base::has_polymermw) {
1955 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1956 const EvalWell water_velocity_abs = abs(water_velocity);
1957 if (poly_inj_conc == 0.) {
1958 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1959 }
1960 const int polymer_table_id = this->polymerTable_();
1961 if (polymer_table_id <= 0) {
1962 OPM_DEFLOG_THROW(std::runtime_error,
1963 fmt::format("Unavailable SKPRPOLY table id used for well {}", name()),
1964 deferred_logger);
1965 }
1966 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1967 const Scalar reference_concentration = skprpolytable.refConcentration;
1968 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1969 // the skin pressure when injecting water, which also means the polymer concentration is zero
1970 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1971 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
1972 if (poly_inj_conc == reference_concentration) {
1973 return sign * pskin_poly;
1974 }
1975 // poly_inj_conc != reference concentration of the table, then some interpolation will be required
1976 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
1977 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
1978 return sign * pskin;
1979 } else {
1980 OPM_DEFLOG_THROW(std::runtime_error,
1981 fmt::format("Polymermw is not activated, while injecting "
1982 "skin pressure is requested for well {}", name()),
1983 deferred_logger);
1984 }
1985 }
1986
1987
1988
1989
1990
1991 template<typename TypeTag>
1994 wpolymermw(const Scalar throughput,
1995 const EvalWell& water_velocity,
1996 DeferredLogger& deferred_logger) const
1997 {
1998 if constexpr (Base::has_polymermw) {
1999 const int table_id = this->polymerInjTable_();
2000 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2001 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2002 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2003 if (this->wpolymer() == 0.) { // not injecting polymer
2004 return molecular_weight;
2005 }
2006 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2007 return molecular_weight;
2008 } else {
2009 OPM_DEFLOG_THROW(std::runtime_error,
2010 fmt::format("Polymermw is not activated, while injecting "
2011 "polymer molecular weight is requested for well {}", name()),
2012 deferred_logger);
2013 }
2014 }
2015
2016
2017
2018
2019
2020 template<typename TypeTag>
2021 void
2023 updateWaterThroughput([[maybe_unused]] const double dt,
2024 WellStateType& well_state) const
2025 {
2026 if constexpr (Base::has_polymermw) {
2027 if (!this->isInjector()) {
2028 return;
2029 }
2030
2031 auto& perf_water_throughput = well_state.well(this->index_of_well_)
2032 .perf_data.water_throughput;
2033
2034 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2035 const Scalar perf_water_vel =
2036 this->primary_variables_.value(Bhp + 1 + perf);
2037
2038 // we do not consider the formation damage due to water
2039 // flowing from reservoir into wellbore
2040 if (perf_water_vel > Scalar{0}) {
2041 perf_water_throughput[perf] += perf_water_vel * dt;
2042 }
2043 }
2044 }
2045 }
2046
2047
2048
2049
2050
2051 template<typename TypeTag>
2052 void
2054 handleInjectivityRate(const Simulator& simulator,
2055 const int perf,
2056 std::vector<EvalWell>& cq_s) const
2057 {
2058 const int cell_idx = this->well_cells_[perf];
2059 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2060 const auto& fs = int_quants.fluidState();
2061 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2062 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2063 const int wat_vel_index = Bhp + 1 + perf;
2064 const unsigned water_comp_idx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx);
2065
2066 // water rate is update to use the form from water velocity, since water velocity is
2067 // a primary variable now
2068 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2069 }
2070
2071
2072
2073
2074 template<typename TypeTag>
2075 void
2077 handleInjectivityEquations(const Simulator& simulator,
2078 const WellStateType& well_state,
2079 const int perf,
2080 const EvalWell& water_flux_s,
2081 DeferredLogger& deferred_logger)
2082 {
2083 const int cell_idx = this->well_cells_[perf];
2084 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2085 const auto& fs = int_quants.fluidState();
2086 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2087 const EvalWell water_flux_r = water_flux_s / b_w;
2088 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2089 const EvalWell water_velocity = water_flux_r / area;
2090 const int wat_vel_index = Bhp + 1 + perf;
2091
2092 // equation for the water velocity
2093 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2094
2095 const auto& ws = well_state.well(this->index_of_well_);
2096 const auto& perf_data = ws.perf_data;
2097 const auto& perf_water_throughput = perf_data.water_throughput;
2098 const Scalar throughput = perf_water_throughput[perf];
2099 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2100
2101 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2102 poly_conc.setValue(this->wpolymer());
2103
2104 // equation for the skin pressure
2105 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2106 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2107
2109 assembleInjectivityEq(eq_pskin,
2110 eq_wat_vel,
2111 pskin_index,
2112 wat_vel_index,
2113 perf,
2114 this->primary_variables_.numWellEq(),
2115 this->linSys_);
2116 }
2117
2118
2119
2120
2121
2122 template<typename TypeTag>
2123 void
2125 checkConvergenceExtraEqs(const std::vector<Scalar>& res,
2126 ConvergenceReport& report) const
2127 {
2128 // if different types of extra equations are involved, this function needs to be refactored further
2129
2130 // checking the convergence of the extra equations related to polymer injectivity
2131 if constexpr (Base::has_polymermw) {
2132 WellConvergence(*this).
2133 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2134 }
2135 }
2136
2137
2138
2139
2140
2141 template<typename TypeTag>
2142 void
2144 updateConnectionRatePolyMW(const EvalWell& cq_s_poly,
2145 const IntensiveQuantities& int_quants,
2146 const WellStateType& well_state,
2147 const int perf,
2148 std::vector<RateVector>& connectionRates,
2149 DeferredLogger& deferred_logger) const
2150 {
2151 // the source term related to transport of molecular weight
2152 EvalWell cq_s_polymw = cq_s_poly;
2153 if (this->isInjector()) {
2154 const int wat_vel_index = Bhp + 1 + perf;
2155 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2156 if (water_velocity > 0.) { // injecting
2157 const auto& ws = well_state.well(this->index_of_well_);
2158 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2159 const Scalar throughput = perf_water_throughput[perf];
2160 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2161 cq_s_polymw *= molecular_weight;
2162 } else {
2163 // we do not consider the molecular weight from the polymer
2164 // going-back to the wellbore through injector
2165 cq_s_polymw *= 0.;
2166 }
2167 } else if (this->isProducer()) {
2168 if (cq_s_polymw < 0.) {
2169 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2170 } else {
2171 // we do not consider the molecular weight from the polymer
2172 // re-injecting back through producer
2173 cq_s_polymw *= 0.;
2174 }
2175 }
2176 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2177 }
2178
2179
2180
2181
2182
2183 template<typename TypeTag>
2184 std::optional<typename StandardWell<TypeTag>::Scalar>
2187 const Simulator& simulator,
2188 const SummaryState& summary_state,
2189 DeferredLogger& deferred_logger) const
2190 {
2191 return computeBhpAtThpLimitProdWithAlq(simulator,
2192 summary_state,
2193 this->getALQ(well_state),
2194 deferred_logger,
2195 /*iterate_if_no_solution */ true);
2196 }
2197
2198 template<typename TypeTag>
2199 std::optional<typename StandardWell<TypeTag>::Scalar>
2202 const SummaryState& summary_state,
2203 const Scalar alq_value,
2204 DeferredLogger& deferred_logger,
2205 bool iterate_if_no_solution) const
2206 {
2207 OPM_TIMEFUNCTION();
2208 // Make the frates() function.
2209 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2210 // Not solving the well equations here, which means we are
2211 // calculating at the current Fg/Fw values of the
2212 // well. This does not matter unless the well is
2213 // crossflowing, and then it is likely still a good
2214 // approximation.
2215 std::vector<Scalar> rates(3);
2216 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2217 this->adaptRatesForVFP(rates);
2218 return rates;
2219 };
2220
2221 Scalar max_pressure = 0.0;
2222 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2223 const int cell_idx = this->well_cells_[perf];
2224 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2225 const auto& fs = int_quants.fluidState();
2226 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2227 max_pressure = std::max(max_pressure, pressure_cell);
2228 }
2229 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2230 summary_state,
2231 max_pressure,
2232 this->connections_.rho(),
2233 alq_value,
2234 this->getTHPConstraint(summary_state),
2235 deferred_logger);
2236
2237 if (bhpAtLimit) {
2238 auto v = frates(*bhpAtLimit);
2239 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2240 return bhpAtLimit;
2241 }
2242 }
2243
2244 if (!iterate_if_no_solution)
2245 return std::nullopt;
2246
2247 auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
2248 // Solver the well iterations to see if we are
2249 // able to get a solution with an update
2250 // solution
2251 std::vector<Scalar> rates(3);
2252 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2253 this->adaptRatesForVFP(rates);
2254 return rates;
2255 };
2256
2257 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2258 summary_state,
2259 max_pressure,
2260 this->connections_.rho(),
2261 alq_value,
2262 this->getTHPConstraint(summary_state),
2263 deferred_logger);
2264
2265
2266 if (bhpAtLimit) {
2267 // should we use fratesIter here since fratesIter is used in computeBhpAtThpLimitProd above?
2268 auto v = frates(*bhpAtLimit);
2269 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2270 return bhpAtLimit;
2271 }
2272 }
2273
2274 // we still don't get a valied solution.
2275 return std::nullopt;
2276 }
2277
2278
2279
2280 template<typename TypeTag>
2281 std::optional<typename StandardWell<TypeTag>::Scalar>
2283 computeBhpAtThpLimitInj(const Simulator& simulator,
2284 const SummaryState& summary_state,
2285 DeferredLogger& deferred_logger) const
2286 {
2287 // Make the frates() function.
2288 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2289 // Not solving the well equations here, which means we are
2290 // calculating at the current Fg/Fw values of the
2291 // well. This does not matter unless the well is
2292 // crossflowing, and then it is likely still a good
2293 // approximation.
2294 std::vector<Scalar> rates(3);
2295 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2296 return rates;
2297 };
2298
2299 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2300 summary_state,
2301 this->connections_.rho(),
2302 1e-6,
2303 50,
2304 true,
2305 deferred_logger);
2306 }
2307
2308
2309
2310
2311
2312 template<typename TypeTag>
2313 bool
2315 iterateWellEqWithControl(const Simulator& simulator,
2316 const double dt,
2317 const Well::InjectionControls& inj_controls,
2318 const Well::ProductionControls& prod_controls,
2319 WellStateType& well_state,
2320 const GroupState<Scalar>& group_state,
2321 DeferredLogger& deferred_logger)
2322 {
2323 updatePrimaryVariables(simulator, well_state, deferred_logger);
2324
2325 const int max_iter = this->param_.max_inner_iter_wells_;
2326 int it = 0;
2327 bool converged;
2328 bool relax_convergence = false;
2329 this->regularize_ = false;
2330 do {
2331 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2332
2333 if (it > this->param_.strict_inner_iter_wells_) {
2334 relax_convergence = true;
2335 this->regularize_ = true;
2336 }
2337
2338 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2339
2340 converged = report.converged();
2341 if (converged) {
2342 break;
2343 }
2344
2345 ++it;
2346 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2347
2348 // TODO: when this function is used for well testing purposes, will need to check the controls, so that we will obtain convergence
2349 // under the most restrictive control. Based on this converged results, we can check whether to re-open the well. Either we refactor
2350 // this function or we use different functions for the well testing purposes.
2351 // We don't allow for switching well controls while computing well potentials and testing wells
2352 // updateWellControl(simulator, well_state, deferred_logger);
2353 } while (it < max_iter);
2354
2355 return converged;
2356 }
2357
2358
2359 template<typename TypeTag>
2360 bool
2362 iterateWellEqWithSwitching(const Simulator& simulator,
2363 const double dt,
2364 const Well::InjectionControls& inj_controls,
2365 const Well::ProductionControls& prod_controls,
2366 WellStateType& well_state,
2367 const GroupState<Scalar>& group_state,
2368 DeferredLogger& deferred_logger,
2369 const bool fixed_control /*false*/,
2370 const bool fixed_status /*false*/)
2371 {
2372 updatePrimaryVariables(simulator, well_state, deferred_logger);
2373
2374 const int max_iter = this->param_.max_inner_iter_wells_;
2375 int it = 0;
2376 bool converged = false;
2377 bool relax_convergence = false;
2378 this->regularize_ = false;
2379 const auto& summary_state = simulator.vanguard().summaryState();
2380
2381 // Always take a few (more than one) iterations after a switch before allowing a new switch
2382 // The optimal number here is subject to further investigation, but it has been observerved
2383 // that unless this number is >1, we may get stuck in a cycle
2384 constexpr int min_its_after_switch = 4;
2385 // We also want to restrict the number of status switches to avoid oscillation between STOP<->OPEN
2386 const int max_status_switch = this->param_.max_well_status_switch_;
2387 int its_since_last_switch = min_its_after_switch;
2388 int switch_count= 0;
2389 // if we fail to solve eqs, we reset status/operability before leaving
2390 const auto well_status_orig = this->wellStatus_;
2391 const auto operability_orig = this->operability_status_;
2392 auto well_status_cur = well_status_orig;
2393 int status_switch_count = 0;
2394 // don't allow opening wells that has a stopped well status
2395 const bool allow_open = well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2396 // don't allow switcing for wells under zero rate target or requested fixed status and control
2397 const bool allow_switching =
2398 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2399 (!fixed_control || !fixed_status) && allow_open;
2400
2401 bool changed = false;
2402 bool final_check = false;
2403 // well needs to be set operable or else solving/updating of re-opened wells is skipped
2404 this->operability_status_.resetOperability();
2405 this->operability_status_.solvable = true;
2406 do {
2407 its_since_last_switch++;
2408 if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
2409 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2410 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2411 inj_controls, prod_controls, wqTotal,
2412 deferred_logger, fixed_control, fixed_status);
2413 if (changed){
2414 its_since_last_switch = 0;
2415 switch_count++;
2416 if (well_status_cur != this->wellStatus_) {
2417 well_status_cur = this->wellStatus_;
2418 status_switch_count++;
2419 }
2420 }
2421 if (!changed && final_check) {
2422 break;
2423 } else {
2424 final_check = false;
2425 }
2426 if (status_switch_count == max_status_switch) {
2427 this->wellStatus_ = well_status_orig;
2428 }
2429 }
2430
2431 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2432
2433 if (it > this->param_.strict_inner_iter_wells_) {
2434 relax_convergence = true;
2435 this->regularize_ = true;
2436 }
2437
2438 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2439
2440 converged = report.converged();
2441 if (converged) {
2442 // if equations are sufficiently linear they might converge in less than min_its_after_switch
2443 // in this case, make sure all constraints are satisfied before returning
2444 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2445 final_check = true;
2446 its_since_last_switch = min_its_after_switch;
2447 } else {
2448 break;
2449 }
2450 }
2451
2452 ++it;
2453 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2454
2455 } while (it < max_iter);
2456
2457 if (converged) {
2458 if (allow_switching){
2459 // update operability if status change
2460 const bool is_stopped = this->wellIsStopped();
2461 if (this->wellHasTHPConstraints(summary_state)){
2462 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2463 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2464 } else {
2465 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2466 }
2467 }
2468 } else {
2469 this->wellStatus_ = well_status_orig;
2470 this->operability_status_ = operability_orig;
2471 const std::string message = fmt::format(" Well {} did not converge in {} inner iterations ("
2472 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2473 deferred_logger.debug(message);
2474 // add operability here as well ?
2475 }
2476 return converged;
2477 }
2478
2479 template<typename TypeTag>
2480 std::vector<typename StandardWell<TypeTag>::Scalar>
2482 computeCurrentWellRates(const Simulator& simulator,
2483 DeferredLogger& deferred_logger) const
2484 {
2485 // Calculate the rates that follow from the current primary variables.
2486 std::vector<Scalar> well_q_s(this->num_conservation_quantities_, 0.);
2487 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2488 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2489 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2490 const int cell_idx = this->well_cells_[perf];
2491 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2492 std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
2493 getMobility(simulator, perf, mob, deferred_logger);
2494 std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
2495 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2496 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2497 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2498 PerforationRates<Scalar> perf_rates;
2499 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2500 cq_s, perf_rates, deferred_logger);
2501 for (int comp = 0; comp < this->num_conservation_quantities_; ++comp) {
2502 well_q_s[comp] += cq_s[comp];
2503 }
2504 }
2505 const auto& comm = this->parallel_well_info_.communication();
2506 if (comm.size() > 1)
2507 {
2508 comm.sum(well_q_s.data(), well_q_s.size());
2509 }
2510 return well_q_s;
2511 }
2512
2513
2514
2515 template <typename TypeTag>
2516 std::vector<typename StandardWell<TypeTag>::Scalar>
2518 getPrimaryVars() const
2519 {
2520 const int num_pri_vars = this->primary_variables_.numWellEq();
2521 std::vector<Scalar> retval(num_pri_vars);
2522 for (int ii = 0; ii < num_pri_vars; ++ii) {
2523 retval[ii] = this->primary_variables_.value(ii);
2524 }
2525 return retval;
2526 }
2527
2528
2529
2530
2531
2532 template <typename TypeTag>
2533 int
2535 setPrimaryVars(typename std::vector<Scalar>::const_iterator it)
2536 {
2537 const int num_pri_vars = this->primary_variables_.numWellEq();
2538 for (int ii = 0; ii < num_pri_vars; ++ii) {
2539 this->primary_variables_.setValue(ii, it[ii]);
2540 }
2541 return num_pri_vars;
2542 }
2543
2544
2545 template <typename TypeTag>
2548 connectionRateEnergy(const std::vector<EvalWell>& cq_s,
2549 const IntensiveQuantities& intQuants,
2550 DeferredLogger& deferred_logger) const
2551 {
2552 auto fs = intQuants.fluidState();
2553 Eval result = 0;
2554 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2555 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2556 continue;
2557 }
2558
2559 // convert to reservoir conditions
2560 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2561 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
2562 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2563 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2564 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2565 } else {
2566 // remove dissolved gas and vapporized oil
2567 const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx);
2568 const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
2569 // q_os = q_or * b_o + rv * q_gr * b_g
2570 // q_gs = q_gr * g_g + rs * q_or * b_o
2571 // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
2572 // d = 1.0 - rs * rv
2573 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2574 if (d <= 0.0) {
2575 deferred_logger.debug(
2576 fmt::format("Problematic d value {} obtained for well {}"
2577 " during calculateSinglePerf with rs {}"
2578 ", rv {}. Continue as if no dissolution (rs = 0) and"
2579 " vaporization (rv = 0) for this connection.",
2580 d, this->name(), fs.Rs(), fs.Rv()));
2581 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2582 } else {
2583 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2584 cq_r_thermal = (cq_s[gasCompIdx] -
2585 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2586 (d * this->extendEval(fs.invB(phaseIdx)) );
2587 } else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2588 // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
2589 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2590 cq_s[gasCompIdx]) /
2591 (d * this->extendEval(fs.invB(phaseIdx)) );
2592 }
2593 }
2594 }
2595
2596 // change temperature for injecting fluids
2597 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2598 // only handles single phase injection now
2599 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2600 fs.setTemperature(this->well_ecl_.inj_temperature());
2601 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
2602 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2603 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2604 paramCache.setRegionIndex(pvtRegionIdx);
2605 paramCache.updatePhase(fs, phaseIdx);
2606
2607 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2608 fs.setDensity(phaseIdx, rho);
2609 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2610 fs.setEnthalpy(phaseIdx, h);
2611 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2612 result += getValue(cq_r_thermal);
2613 } else if (cq_r_thermal > 0.0) {
2614 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2615 result += Base::restrictEval(cq_r_thermal);
2616 } else {
2617 // compute the thermal flux
2618 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2619 result += Base::restrictEval(cq_r_thermal);
2620 }
2621 }
2622
2623 return result * this->well_efficiency_factor_;
2624 }
2625} // namespace Opm
2626
2627#endif
#define OPM_DEFLOG_THROW(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:45
#define OPM_DEFLOG_PROBLEM(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:61
Definition: ConvergenceReport.hpp:38
Definition: DeferredLogger.hpp:57
void warning(const std::string &tag, const std::string &message)
void debug(const std::string &tag, const std::string &message)
Definition: GroupState.hpp:41
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:198
Definition: RatioCalculator.hpp:38
Class handling assemble of the equation system for StandardWell.
Definition: StandardWellAssemble.hpp:43
Scalar pressure_diff(const unsigned perf) const
Returns pressure drop for a given perforation.
Definition: StandardWellConnections.hpp:106
StdWellConnections connections_
Connection level values.
Definition: StandardWellEval.hpp:105
PrimaryVariables primary_variables_
Primary variables for well.
Definition: StandardWellEval.hpp:99
Definition: StandardWell.hpp:60
void calculateExplicitQuantities(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1380
EvalWell wpolymermw(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1994
typename StdWellEval::EvalWell EvalWell
Definition: StandardWell.hpp:119
void updateWellStateFromPrimaryVariables(WellStateType &well_state, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:760
WellConnectionProps computePropertiesForWellConnectionPressures(const Simulator &simulator, const WellStateType &well_state) const
Definition: StandardWell_impl.hpp:1119
typename StdWellEval::BVectorWell BVectorWell
Definition: StandardWell.hpp:120
std::optional< Scalar > computeBhpAtThpLimitProd(const WellStateType &well_state, const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2186
void computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1699
void addWellContributions(SparseMatrixAdapter &mat) const override
Definition: StandardWell_impl.hpp:1889
std::vector< Scalar > getPrimaryVars() const override
Definition: StandardWell_impl.hpp:2518
void addWellPressureEquations(PressureMatrix &mat, const BVector &x, const int pressureVarIndex, const bool use_well_weights, const WellStateType &well_state) const override
Definition: StandardWell_impl.hpp:1897
void updateWaterMobilityWithPolymer(const Simulator &simulator, const int perf, std::vector< EvalWell > &mob_water, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1825
void assembleWellEqWithoutIteration(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:341
std::vector< Scalar > computeCurrentWellRates(const Simulator &ebosSimulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:2482
void calculateSinglePerf(const Simulator &simulator, const int perf, WellStateType &well_state, std::vector< RateVector > &connectionRates, std::vector< EvalWell > &cq_s, EvalWell &water_flux_s, EvalWell &cq_s_zfrac_effective, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:493
void updatePrimaryVariablesNewton(const BVectorWell &dwells, const bool stop_or_zero_rate_target, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:737
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:82
void computeWellConnectionPressures(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1340
StandardWell(const Well &well, const ParallelWellInfo< Scalar > &pw_info, const int time_step, const ModelParameters &param, const RateConverterType &rate_converter, const int pvtRegionIdx, const int num_conservation_quantities, const int num_phases, const int index_of_well, const std::vector< PerforationData< Scalar > > &perf_data)
Definition: StandardWell_impl.hpp:52
typename StdWellEval::StdWellConnections::Properties WellConnectionProps
Definition: StandardWell.hpp:263
void updateConnectionRatePolyMW(const EvalWell &cq_s_poly, const IntensiveQuantities &int_quants, const WellStateType &well_state, const int perf, std::vector< RateVector > &connectionRates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2144
bool iterateWellEqWithSwitching(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger, const bool fixed_control=false, const bool fixed_status=false) override
Definition: StandardWell_impl.hpp:2362
void computeWellRatesWithBhp(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1446
void getMobility(const Simulator &simulator, const int perf, std::vector< Value > &mob, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:659
void computeWellRatesWithBhpIterations(const Simulator &ebosSimulator, const Scalar &bhp, std::vector< Scalar > &well_flux, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1492
void updateIPR(const Simulator &simulator, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:779
std::optional< Scalar > computeBhpAtThpLimitProdWithAlq(const Simulator &ebos_simulator, const SummaryState &summary_state, const Scalar alq_value, DeferredLogger &deferred_logger, bool iterate_if_no_solution) const override
Definition: StandardWell_impl.hpp:2201
void updateIPRImplicit(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:872
void computeWellConnectionDensitesPressures(const Simulator &simulator, const WellStateType &well_state, const WellConnectionProps &props, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:1289
void computePerfRate(const IntensiveQuantities &intQuants, const std::vector< Value > &mob, const Value &bhp, const std::vector< Scalar > &Tw, const int perf, const bool allow_cf, std::vector< Value > &cq_s, PerforationRates< Scalar > &perf_rates, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:93
void updateWellState(const Simulator &simulator, const BVectorWell &dwells, WellStateType &well_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:715
void handleInjectivityEquations(const Simulator &simulator, const WellStateType &well_state, const int perf, const EvalWell &water_flux_s, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:2077
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1394
virtual ConvergenceReport getWellConvergence(const Simulator &simulator, const WellStateType &well_state, const std::vector< Scalar > &B_avg, DeferredLogger &deferred_logger, const bool relax_tolerance) const override
check whether the well equations get converged for this well
Definition: StandardWell_impl.hpp:1170
void checkConvergenceExtraEqs(const std::vector< Scalar > &res, ConvergenceReport &report) const
Definition: StandardWell_impl.hpp:2125
typename StdWellEval::Eval Eval
Definition: StandardWell.hpp:118
Scalar computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, std::vector< Scalar > &potentials, Scalar alq) const
Definition: StandardWell_impl.hpp:1670
bool openCrossFlowAvoidSingularity(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1108
std::optional< Scalar > computeBhpAtThpLimitInj(const Simulator &simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:2283
bool allDrawDownWrongDirection(const Simulator &simulator) const
Definition: StandardWell_impl.hpp:1066
EvalWell pskin(const Scalar throughput, const EvalWell &water_velocity, const EvalWell &poly_inj_conc, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1949
static constexpr int numWellConservationEq
Definition: StandardWell.hpp:95
int setPrimaryVars(typename std::vector< Scalar >::const_iterator it) override
Definition: StandardWell_impl.hpp:2535
void checkOperabilityUnderTHPLimit(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1017
void computeWellPotentials(const Simulator &simulator, const WellStateType &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) override
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1716
bool computeWellPotentialsImplicit(const Simulator &ebos_simulator, const WellStateType &well_state, std::vector< Scalar > &well_potentials, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1591
void updateWaterThroughput(const double dt, WellStateType &well_state) const override
Definition: StandardWell_impl.hpp:2023
void checkOperabilityUnderBHPLimit(const WellStateType &well_state, const Simulator &simulator, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:948
void recoverWellSolutionAndUpdateWellState(const Simulator &simulator, const BVector &x, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1426
void assembleWellEqWithoutIterationImpl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: StandardWell_impl.hpp:367
bool iterateWellEqWithControl(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:2315
EvalWell pskinwater(const Scalar throughput, const EvalWell &water_velocity, DeferredLogger &deferred_logger) const
Definition: StandardWell_impl.hpp:1917
void solveEqAndUpdateWellState(const Simulator &simulator, WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1358
void handleInjectivityRate(const Simulator &simulator, const int perf, std::vector< EvalWell > &cq_s) const
Definition: StandardWell_impl.hpp:2054
virtual void init(const std::vector< Scalar > &depth_arg, const Scalar gravity_arg, const std::vector< Scalar > &B_avg, const bool changed_to_open_this_step) override
Definition: StandardWell_impl.hpp:76
std::vector< Scalar > computeWellPotentialWithTHP(const Simulator &ebosSimulator, DeferredLogger &deferred_logger, const WellStateType &well_state) const
Definition: StandardWell_impl.hpp:1556
void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator< Scalar > &wellPICalc, WellStateType &well_state, DeferredLogger &deferred_logger) const override
Definition: StandardWell_impl.hpp:1214
void updatePrimaryVariables(const Simulator &simulator, const WellStateType &well_state, DeferredLogger &deferred_logger) override
Definition: StandardWell_impl.hpp:1791
Scalar getRefDensity() const override
Definition: StandardWell_impl.hpp:1814
Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const override
Definition: StandardWell_impl.hpp:1776
EvalWell getQs(const int compIdx) const
Returns scaled rate for a component.
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
Scalar calculateThpFromBhp(const std::vector< Scalar > &rates, const Scalar bhp, const Scalar rho, const std::optional< Scalar > &alq, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Calculates THP from BHP.
std::optional< Scalar > computeBhpAtThpLimitProd(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar maxPerfPress, const Scalar rho, const Scalar alq_value, const Scalar thp_limit, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for a producer.
Scalar mostStrictBhpFromBhpLimits(const SummaryState &summaryState) const
Obtain the most strict BHP from BHP limits.
std::optional< Scalar > computeBhpAtThpLimitInj(const std::function< std::vector< Scalar >(const Scalar)> &frates, const SummaryState &summary_state, const Scalar rho, const Scalar flo_rel_tol, const int max_iteration, const bool throwOnError, DeferredLogger &deferred_logger) const
Compute BHP from THP limit for an injector.
Definition: WellConvergence.hpp:38
const int num_conservation_quantities_
Definition: WellInterfaceGeneric.hpp:310
Well well_ecl_
Definition: WellInterfaceGeneric.hpp:300
void resetDampening()
Definition: WellInterfaceGeneric.hpp:243
std::pair< bool, bool > computeWellPotentials(std::vector< Scalar > &well_potentials, const WellStateType &well_state)
void prepareForPotentialCalculations(const SummaryState &summary_state, WellStateType &well_state, Well::InjectionControls &inj_controls, Well::ProductionControls &prod_controls) const
Definition: WellInterfaceIndices.hpp:34
Definition: WellInterface.hpp:76
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: WellInterface.hpp:81
typename WellInterfaceFluidSystem< FluidSystem >::RateConverterType RateConverterType
Definition: WellInterface.hpp:103
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, 1, 1 > > PressureMatrix
Definition: WellInterface.hpp:97
void getMobility(const Simulator &simulator, const int local_perf_index, std::vector< Value > &mob, Callback &extendEval, DeferredLogger &deferred_logger) const
Definition: WellInterface_impl.hpp:1899
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: WellInterface.hpp:86
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: WellInterface.hpp:82
Dune::BlockVector< VectorBlockType > BVector
Definition: WellInterface.hpp:96
typename Base::ModelParameters ModelParameters
Definition: WellInterface.hpp:109
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: WellInterface.hpp:83
bool solveWellWithOperabilityCheck(const Simulator &simulator, const double dt, const Well::InjectionControls &inj_controls, const Well::ProductionControls &prod_controls, WellStateType &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:563
GetPropType< TypeTag, Properties::Indices > Indices
Definition: WellInterface.hpp:85
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: WellInterface.hpp:88
Definition: WellProdIndexCalculator.hpp:37
Scalar connectionProdIndStandard(const std::size_t connIdx, const Scalar connMobility) const
Definition: WellState.hpp:66
const SingleWellState< Scalar, IndexTraits > & well(std::size_t well_index) const
Definition: WellState.hpp:289
std::vector< Scalar > & wellRates(std::size_t well_index)
One rate per well and phase.
Definition: WellState.hpp:254
@ NONE
Definition: DeferredLogger.hpp:46
Definition: blackoilboundaryratevector.hh:39
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Static data associated with a well perforation.
Definition: PerforationData.hpp:30
Definition: PerforationData.hpp:41
Scalar dis_gas
Definition: PerforationData.hpp:42
Scalar vap_wat
Definition: PerforationData.hpp:45
Scalar vap_oil
Definition: PerforationData.hpp:44
Scalar dis_gas_in_water
Definition: PerforationData.hpp:43