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