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