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