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