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