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