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