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(perf_data.microbial_rates[perf],
627 perf_data.oxygen_rates[perf],
628 perf_data.urea_rates[perf],
629 cq_s,
630 microbialConcentration,
631 oxygenConcentration,
632 ureaConcentration);
633 }
634
635 // Store the perforation pressure for later usage.
636 perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
637
638 // Store the perforation gass mass rate.
639 const auto& pu = well_state.phaseUsage();
640 if (pu.has_co2_or_h2store) {
641 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
642 const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
643 perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
644 }
645 }
646
647
648
649 template<typename TypeTag>
650 template<class Value>
651 void
652 StandardWell<TypeTag>::
653 getMobility(const Simulator& simulator,
654 const int perf,
655 std::vector<Value>& mob,
656 DeferredLogger& deferred_logger) const
657 {
658 auto obtain = [this](const Eval& value)
659 {
660 if constexpr (std::is_same_v<Value, Scalar>) {
661 static_cast<void>(this); // suppress clang warning
662 return getValue(value);
663 } else {
664 return this->extendEval(value);
665 }
666 };
667 WellInterface<TypeTag>::getMobility(simulator, perf, mob,
668 obtain, deferred_logger);
669
670 // modify the water mobility if polymer is present
671 if constexpr (has_polymer) {
672 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
673 OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
674 }
675
676 // for the cases related to polymer molecular weight, we assume fully mixing
677 // as a result, the polymer and water share the same viscosity
678 if constexpr (!Base::has_polymermw) {
679 if constexpr (std::is_same_v<Value, Scalar>) {
680 std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
681 for (std::size_t i = 0; i < mob.size(); ++i) {
682 mob_eval[i].setValue(mob[i]);
683 }
684 updateWaterMobilityWithPolymer(simulator, perf, mob_eval, deferred_logger);
685 for (std::size_t i = 0; i < mob.size(); ++i) {
686 mob[i] = getValue(mob_eval[i]);
687 }
688 } else {
689 updateWaterMobilityWithPolymer(simulator, perf, mob, deferred_logger);
690 }
691 }
692 }
693
694 // if the injecting well has WINJMULT setup, we update the mobility accordingly
695 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
696 const Scalar bhp = this->primary_variables_.value(Bhp);
697 const Scalar perf_press = bhp + this->connections_.pressure_diff(perf);
698 const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
699 for (std::size_t i = 0; i < mob.size(); ++i) {
700 mob[i] *= multiplier;
701 }
702 }
703 }
704
705
706 template<typename TypeTag>
707 void
708 StandardWell<TypeTag>::
709 updateWellState(const Simulator& simulator,
710 const BVectorWell& dwells,
711 WellState<Scalar>& well_state,
712 DeferredLogger& deferred_logger)
713 {
714 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
715
716 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
717 updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
718
719 const auto& summary_state = simulator.vanguard().summaryState();
720 updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
721 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
722 }
723
724
725
726
727
728 template<typename TypeTag>
729 void
730 StandardWell<TypeTag>::
731 updatePrimaryVariablesNewton(const BVectorWell& dwells,
732 const bool stop_or_zero_rate_target,
733 DeferredLogger& deferred_logger)
734 {
735 const Scalar dFLimit = this->param_.dwell_fraction_max_;
736 const Scalar dBHPLimit = this->param_.dbhp_max_rel_;
737 this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger);
738
739 // for the water velocity and skin pressure
740 if constexpr (Base::has_polymermw) {
741 this->primary_variables_.updateNewtonPolyMW(dwells);
742 }
743
744 this->primary_variables_.checkFinite(deferred_logger);
745 }
746
747
748
749
750
751 template<typename TypeTag>
752 void
753 StandardWell<TypeTag>::
754 updateWellStateFromPrimaryVariables(WellState<Scalar>& well_state,
755 const SummaryState& summary_state,
756 DeferredLogger& deferred_logger) const
757 {
758 this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, summary_state, deferred_logger);
759
760 // other primary variables related to polymer injectivity study
761 if constexpr (Base::has_polymermw) {
762 this->primary_variables_.copyToWellStatePolyMW(well_state);
763 }
764 }
765
766
767
768
769
770 template<typename TypeTag>
771 void
772 StandardWell<TypeTag>::
773 updateIPR(const Simulator& simulator, DeferredLogger& deferred_logger) const
774 {
775 // TODO: not handling solvent related here for now
776
777 // initialize all the values to be zero to begin with
778 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
779 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
780
781 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
782 std::vector<Scalar> mob(this->num_components_, 0.0);
783 getMobility(simulator, perf, mob, deferred_logger);
784
785 const int cell_idx = this->well_cells_[perf];
786 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
787 const auto& fs = int_quantities.fluidState();
788 // the pressure of the reservoir grid block the well connection is in
789 Scalar p_r = this->getPerfCellPressure(fs).value();
790
791 // calculating the b for the connection
792 std::vector<Scalar> b_perf(this->num_components_);
793 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
794 if (!FluidSystem::phaseIsActive(phase)) {
795 continue;
796 }
797 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
798 b_perf[comp_idx] = fs.invB(phase).value();
799 }
800 if constexpr (has_solvent) {
801 b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
802 }
803
804 // the pressure difference between the connection and BHP
805 const Scalar h_perf = this->connections_.pressure_diff(perf);
806 const Scalar pressure_diff = p_r - h_perf;
807
808 // Let us add a check, since the pressure is calculated based on zero value BHP
809 // it should not be negative anyway. If it is negative, we might need to re-formulate
810 // to taking into consideration the crossflow here.
811 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
812 deferred_logger.debug("CROSSFLOW_IPR",
813 "cross flow found when updateIPR for well " + name()
814 + " . The connection is ignored in IPR calculations");
815 // we ignore these connections for now
816 continue;
817 }
818
819 // the well index associated with the connection
820 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
821 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
822 const std::vector<Scalar> tw_perf = this->wellIndex(perf,
823 int_quantities,
824 trans_mult,
825 wellstate_nupcol);
826 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
827 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
828 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
829 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
830 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
831 ipr_b_perf[comp_idx] += tw_mob;
832 }
833
834 // we need to handle the rs and rv when both oil and gas are present
835 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
836 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
837 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
838 const Scalar rs = (fs.Rs()).value();
839 const Scalar rv = (fs.Rv()).value();
840
841 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
842 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
843
844 ipr_a_perf[gas_comp_idx] += dis_gas_a;
845 ipr_a_perf[oil_comp_idx] += vap_oil_a;
846
847 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
848 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
849
850 ipr_b_perf[gas_comp_idx] += dis_gas_b;
851 ipr_b_perf[oil_comp_idx] += vap_oil_b;
852 }
853
854 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
855 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
856 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
857 }
858 }
859 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
860 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
861 }
862
863 template<typename TypeTag>
864 void
865 StandardWell<TypeTag>::
866 updateIPRImplicit(const Simulator& simulator,
867 WellState<Scalar>& well_state,
868 DeferredLogger& deferred_logger)
869 {
870 // Compute IPR based on *converged* well-equation:
871 // For a component rate r the derivative dr/dbhp is obtained by
872 // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target)
873 // where Eq(x)=0 is the well equation setup with bhp control and primary variables x
874
875 // We shouldn't have zero rates at this stage, but check
876 bool zero_rates;
877 auto rates = well_state.well(this->index_of_well_).surface_rates;
878 zero_rates = true;
879 for (std::size_t p = 0; p < rates.size(); ++p) {
880 zero_rates &= rates[p] == 0.0;
881 }
882 auto& ws = well_state.well(this->index_of_well_);
883 if (zero_rates) {
884 const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
885 deferred_logger.debug(msg);
886 /*
887 // could revert to standard approach here:
888 updateIPR(simulator, deferred_logger);
889 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
890 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
891 ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
892 ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
893 }
894 return;
895 */
896 }
897 const auto& group_state = simulator.problem().wellModel().groupState();
898
899 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
900 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
901
902 auto inj_controls = Well::InjectionControls(0);
903 auto prod_controls = Well::ProductionControls(0);
904 prod_controls.addControl(Well::ProducerCMode::BHP);
905 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
906
907 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
908 const auto cmode = ws.production_cmode;
909 ws.production_cmode = Well::ProducerCMode::BHP;
910 const double dt = simulator.timeStepSize();
911 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
912
913 const size_t nEq = this->primary_variables_.numWellEq();
914 BVectorWell rhs(1);
915 rhs[0].resize(nEq);
916 // rhs = 0 except -1 for control eq
917 for (size_t i=0; i < nEq; ++i){
918 rhs[0][i] = 0.0;
919 }
920 rhs[0][Bhp] = -1.0;
921
922 BVectorWell x_well(1);
923 x_well[0].resize(nEq);
924 this->linSys_.solve(rhs, x_well);
925
926 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
927 EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
928 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
929 for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
930 // well primary variable derivatives in EvalWell start at position Indices::numEq
931 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
932 }
933 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
934 }
935 // reset cmode
936 ws.production_cmode = cmode;
937 }
938
939 template<typename TypeTag>
940 void
941 StandardWell<TypeTag>::
942 checkOperabilityUnderBHPLimit(const WellState<Scalar>& well_state,
943 const Simulator& simulator,
944 DeferredLogger& deferred_logger)
945 {
946 const auto& summaryState = simulator.vanguard().summaryState();
947 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
948 // Crude but works: default is one atmosphere.
949 // TODO: a better way to detect whether the BHP is defaulted or not
950 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
951 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
952 // if the BHP limit is not defaulted or the well does not have a THP limit
953 // we need to check the BHP limit
954 Scalar total_ipr_mass_rate = 0.0;
955 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
956 {
957 if (!FluidSystem::phaseIsActive(phaseIdx)) {
958 continue;
959 }
960
961 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
962 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
963
964 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
965 total_ipr_mass_rate += ipr_rate * rho;
966 }
967 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
968 this->operability_status_.operable_under_only_bhp_limit = false;
969 }
970
971 // checking whether running under BHP limit will violate THP limit
972 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
973 // option 1: calculate well rates based on the BHP limit.
974 // option 2: stick with the above IPR curve
975 // we use IPR here
976 std::vector<Scalar> well_rates_bhp_limit;
977 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
978
979 this->adaptRatesForVFP(well_rates_bhp_limit);
980 const Scalar thp_limit = this->getTHPConstraint(summaryState);
981 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
982 bhp_limit,
983 this->connections_.rho(),
984 this->getALQ(well_state),
985 thp_limit,
986 deferred_logger);
987 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
988 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
989 }
990 }
991 } else {
992 // defaulted BHP and there is a THP constraint
993 // default BHP limit is about 1 atm.
994 // when applied the hydrostatic pressure correction dp,
995 // most likely we get a negative value (bhp + dp)to search in the VFP table,
996 // which is not desirable.
997 // we assume we can operate under defaulted BHP limit and will violate the THP limit
998 // when operating under defaulted BHP limit.
999 this->operability_status_.operable_under_only_bhp_limit = true;
1000 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1001 }
1002 }
1003
1004
1005
1006
1007
1008 template<typename TypeTag>
1009 void
1010 StandardWell<TypeTag>::
1011 checkOperabilityUnderTHPLimit(const Simulator& simulator,
1012 const WellState<Scalar>& well_state,
1013 DeferredLogger& deferred_logger)
1014 {
1015 const auto& summaryState = simulator.vanguard().summaryState();
1016 const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger)
1017 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1018
1019 if (obtain_bhp) {
1020 this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1021
1022 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1023 this->operability_status_.obey_bhp_limit_with_thp_limit = this->isProducer() ?
1024 *obtain_bhp >= bhp_limit : *obtain_bhp <= bhp_limit ;
1025
1026 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1027 if (this->isProducer() && *obtain_bhp < thp_limit) {
1028 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1029 + " bars is SMALLER than thp limit "
1030 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1031 + " bars as a producer for well " + name();
1032 deferred_logger.debug(msg);
1033 }
1034 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1035 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1036 + " bars is LARGER than thp limit "
1037 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1038 + " bars as a injector for well " + name();
1039 deferred_logger.debug(msg);
1040 }
1041 } else {
1042 this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1043 this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1044 if (!this->wellIsStopped()) {
1045 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1046 deferred_logger.debug(" could not find bhp value at thp limit "
1047 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1048 + " bar for well " + name() + ", the well might need to be closed ");
1049 }
1050 }
1051 }
1052
1053
1054
1055
1056
1057 template<typename TypeTag>
1058 bool
1059 StandardWell<TypeTag>::
1060 allDrawDownWrongDirection(const Simulator& simulator) const
1061 {
1062 bool all_drawdown_wrong_direction = true;
1063
1064 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1065 const int cell_idx = this->well_cells_[perf];
1066 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1067 const auto& fs = intQuants.fluidState();
1068
1069 const Scalar pressure = this->getPerfCellPressure(fs).value();
1070 const Scalar bhp = this->primary_variables_.eval(Bhp).value();
1071
1072 // Pressure drawdown (also used to determine direction of flow)
1073 const Scalar well_pressure = bhp + this->connections_.pressure_diff(perf);
1074 const Scalar drawdown = pressure - well_pressure;
1075
1076 // for now, if there is one perforation can produce/inject in the correct
1077 // direction, we consider this well can still produce/inject.
1078 // TODO: it can be more complicated than this to cause wrong-signed rates
1079 if ( (drawdown < 0. && this->isInjector()) ||
1080 (drawdown > 0. && this->isProducer()) ) {
1081 all_drawdown_wrong_direction = false;
1082 break;
1083 }
1084 }
1085
1086 const auto& comm = this->parallel_well_info_.communication();
1087 if (comm.size() > 1)
1088 {
1089 all_drawdown_wrong_direction =
1090 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1091 }
1092
1093 return all_drawdown_wrong_direction;
1094 }
1095
1096
1097
1098
1099 template<typename TypeTag>
1100 bool
1101 StandardWell<TypeTag>::
1102 canProduceInjectWithCurrentBhp(const Simulator& simulator,
1103 const WellState<Scalar>& well_state,
1104 DeferredLogger& deferred_logger)
1105 {
1106 const Scalar bhp = well_state.well(this->index_of_well_).bhp;
1107 std::vector<Scalar> well_rates;
1108 computeWellRatesWithBhp(simulator, bhp, well_rates, deferred_logger);
1109
1110 const Scalar sign = (this->isProducer()) ? -1. : 1.;
1111 const Scalar threshold = sign * std::numeric_limits<Scalar>::min();
1112
1113 bool can_produce_inject = false;
1114 for (const auto value : well_rates) {
1115 if (this->isProducer() && value < threshold) {
1116 can_produce_inject = true;
1117 break;
1118 } else if (this->isInjector() && value > threshold) {
1119 can_produce_inject = true;
1120 break;
1121 }
1122 }
1123
1124 if (!can_produce_inject) {
1125 deferred_logger.debug(" well " + name() + " CANNOT produce or inejct ");
1126 }
1127
1128 return can_produce_inject;
1129 }
1130
1131
1132
1133
1134
1135 template<typename TypeTag>
1136 bool
1137 StandardWell<TypeTag>::
1138 openCrossFlowAvoidSingularity(const Simulator& simulator) const
1139 {
1140 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1141 }
1142
1143
1144
1145
1146 template<typename TypeTag>
1147 typename StandardWell<TypeTag>::WellConnectionProps
1148 StandardWell<TypeTag>::
1149 computePropertiesForWellConnectionPressures(const Simulator& simulator,
1150 const WellState<Scalar>& well_state) const
1151 {
1152 auto prop_func = typename StdWellEval::StdWellConnections::PressurePropertyFunctions {
1153 // getTemperature
1154 [&model = simulator.model()](int cell_idx, int phase_idx)
1155 {
1156 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1157 .fluidState().temperature(phase_idx).value();
1158 },
1159
1160 // getSaltConcentration
1161 [&model = simulator.model()](int cell_idx)
1162 {
1163 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1164 .fluidState().saltConcentration().value();
1165 },
1166
1167 // getPvtRegionIdx
1168 [&model = simulator.model()](int cell_idx)
1169 {
1170 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1171 .fluidState().pvtRegionIndex();
1172 }
1173 };
1174
1175 if constexpr (Indices::enableSolvent) {
1176 prop_func.solventInverseFormationVolumeFactor =
1177 [&model = simulator.model()](int cell_idx)
1178 {
1179 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1180 .solventInverseFormationVolumeFactor().value();
1181 };
1182
1183 prop_func.solventRefDensity = [&model = simulator.model()](int cell_idx)
1184 {
1185 return model.intensiveQuantities(cell_idx, /* time_idx = */ 0)
1186 .solventRefDensity();
1187 };
1188 }
1189
1190 return this->connections_.computePropertiesForPressures(well_state, prop_func);
1191 }
1192
1193
1194
1195
1196
1197 template<typename TypeTag>
1198 ConvergenceReport
1200 getWellConvergence(const Simulator& simulator,
1201 const WellState<Scalar>& well_state,
1202 const std::vector<Scalar>& B_avg,
1203 DeferredLogger& deferred_logger,
1204 const bool relax_tolerance) const
1205 {
1206 // the following implementation assume that the polymer is always after the w-o-g phases
1207 // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
1208 assert((int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1209
1210 Scalar tol_wells = this->param_.tolerance_wells_;
1211 // use stricter tolerance for stopped wells and wells under zero rate target control.
1212 constexpr Scalar stopped_factor = 1.e-4;
1213 // use stricter tolerance for dynamic thp to ameliorate network convergence
1214 constexpr Scalar dynamic_thp_factor = 1.e-1;
1215 if (this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger)) {
1216 tol_wells = tol_wells*stopped_factor;
1217 } else if (this->getDynamicThpLimit()) {
1218 tol_wells = tol_wells*dynamic_thp_factor;
1219 }
1220
1221 std::vector<Scalar> res;
1222 ConvergenceReport report = this->StdWellEval::getWellConvergence(well_state,
1223 B_avg,
1224 this->param_.max_residual_allowed_,
1225 tol_wells,
1226 this->param_.relaxed_tolerance_flow_well_,
1227 relax_tolerance,
1228 this->wellIsStopped(),
1229 res,
1230 deferred_logger);
1231
1232 checkConvergenceExtraEqs(res, report);
1233
1234 return report;
1235 }
1236
1237
1238
1239
1240
1241 template<typename TypeTag>
1242 void
1244 updateProductivityIndex(const Simulator& simulator,
1245 const WellProdIndexCalculator<Scalar>& wellPICalc,
1246 WellState<Scalar>& well_state,
1247 DeferredLogger& deferred_logger) const
1248 {
1249 auto fluidState = [&simulator, this](const int perf)
1250 {
1251 const auto cell_idx = this->well_cells_[perf];
1252 return simulator.model()
1253 .intensiveQuantities(cell_idx, /*timeIdx=*/ 0).fluidState();
1254 };
1255
1256 const int np = this->number_of_phases_;
1257 auto setToZero = [np](Scalar* x) -> void
1258 {
1259 std::fill_n(x, np, 0.0);
1260 };
1261
1262 auto addVector = [np](const Scalar* src, Scalar* dest) -> void
1263 {
1264 std::transform(src, src + np, dest, dest, std::plus<>{});
1265 };
1266
1267 auto& ws = well_state.well(this->index_of_well_);
1268 auto& perf_data = ws.perf_data;
1269 auto* wellPI = ws.productivity_index.data();
1270 auto* connPI = perf_data.prod_index.data();
1271
1272 setToZero(wellPI);
1273
1274 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1275 auto subsetPerfID = 0;
1276
1277 for (const auto& perf : *this->perf_data_) {
1278 auto allPerfID = perf.ecl_index;
1279
1280 auto connPICalc = [&wellPICalc, allPerfID](const Scalar mobility) -> Scalar
1281 {
1282 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
1283 };
1284
1285 std::vector<Scalar> mob(this->num_components_, 0.0);
1286 getMobility(simulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
1287
1288 const auto& fs = fluidState(subsetPerfID);
1289 setToZero(connPI);
1290
1291 if (this->isInjector()) {
1292 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1293 mob, connPI, deferred_logger);
1294 }
1295 else { // Production or zero flow rate
1296 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1297 }
1298
1299 addVector(connPI, wellPI);
1300
1301 ++subsetPerfID;
1302 connPI += np;
1303 }
1304
1305 // Sum with communication in case of distributed well.
1306 const auto& comm = this->parallel_well_info_.communication();
1307 if (comm.size() > 1) {
1308 comm.sum(wellPI, np);
1309 }
1310
1311 assert ((static_cast<int>(subsetPerfID) == this->number_of_local_perforations_) &&
1312 "Internal logic error in processing connections for PI/II");
1313 }
1314
1315
1316
1317 template<typename TypeTag>
1318 void StandardWell<TypeTag>::
1319 computeWellConnectionDensitesPressures(const Simulator& simulator,
1320 const WellState<Scalar>& well_state,
1321 const WellConnectionProps& props,
1322 DeferredLogger& deferred_logger)
1323 {
1324 // Cell level dynamic property call-back functions as fall-back
1325 // option for calculating connection level mixture densities in
1326 // stopped or zero-rate producer wells.
1327 const auto prop_func = typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
1328 // This becomes slightly more palatable with C++20's designated
1329 // initialisers.
1330
1331 // mobility: Phase mobilities in specified cell.
1332 [&model = simulator.model()](const int cell,
1333 const std::vector<int>& phases,
1334 std::vector<Scalar>& mob)
1335 {
1336 const auto& iq = model.intensiveQuantities(cell, /* time_idx = */ 0);
1337
1338 std::transform(phases.begin(), phases.end(), mob.begin(),
1339 [&iq](const int phase) { return iq.mobility(phase).value(); });
1340 },
1341
1342 // densityInCell: Reservoir condition phase densities in
1343 // specified cell.
1344 [&model = simulator.model()](const int cell,
1345 const std::vector<int>& phases,
1346 std::vector<Scalar>& rho)
1347 {
1348 const auto& fs = model.intensiveQuantities(cell, /* time_idx = */ 0).fluidState();
1349
1350 std::transform(phases.begin(), phases.end(), rho.begin(),
1351 [&fs](const int phase) { return fs.density(phase).value(); });
1352 }
1353 };
1354
1355 const auto stopped_or_zero_rate_target = this->
1356 stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1357
1358 this->connections_
1359 .computeProperties(stopped_or_zero_rate_target, well_state,
1360 prop_func, props, deferred_logger);
1361 }
1362
1363
1364
1365
1366
1367 template<typename TypeTag>
1368 void
1369 StandardWell<TypeTag>::
1370 computeWellConnectionPressures(const Simulator& simulator,
1371 const WellState<Scalar>& well_state,
1372 DeferredLogger& deferred_logger)
1373 {
1374 const auto props = computePropertiesForWellConnectionPressures
1375 (simulator, well_state);
1376
1377 computeWellConnectionDensitesPressures(simulator, well_state,
1378 props, deferred_logger);
1379 }
1380
1381
1382
1383
1384
1385 template<typename TypeTag>
1386 void
1387 StandardWell<TypeTag>::
1388 solveEqAndUpdateWellState(const Simulator& simulator,
1389 WellState<Scalar>& well_state,
1390 DeferredLogger& deferred_logger)
1391 {
1392 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1393
1394 // We assemble the well equations, then we check the convergence,
1395 // which is why we do not put the assembleWellEq here.
1396 BVectorWell dx_well(1);
1397 dx_well[0].resize(this->primary_variables_.numWellEq());
1398 this->linSys_.solve( dx_well);
1399
1400 updateWellState(simulator, dx_well, well_state, deferred_logger);
1401 }
1402
1403
1404
1405
1406
1407 template<typename TypeTag>
1408 void
1409 StandardWell<TypeTag>::
1410 calculateExplicitQuantities(const Simulator& simulator,
1411 const WellState<Scalar>& well_state,
1412 DeferredLogger& deferred_logger)
1413 {
1414 updatePrimaryVariables(simulator, well_state, deferred_logger);
1415 computeWellConnectionPressures(simulator, well_state, deferred_logger);
1416 this->computeAccumWell();
1417 }
1418
1419
1420
1421 template<typename TypeTag>
1422 void
1424 apply(const BVector& x, BVector& Ax) const
1425 {
1426 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1427
1428 if (this->param_.matrix_add_well_contributions_)
1429 {
1430 // Contributions are already in the matrix itself
1431 return;
1432 }
1433
1434 this->linSys_.apply(x, Ax);
1435 }
1436
1437
1438
1439
1440 template<typename TypeTag>
1441 void
1443 apply(BVector& r) const
1444 {
1445 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1446
1447 this->linSys_.apply(r);
1448 }
1449
1450
1451
1452
1453 template<typename TypeTag>
1454 void
1456 recoverWellSolutionAndUpdateWellState(const Simulator& simulator,
1457 const BVector& x,
1458 WellState<Scalar>& well_state,
1459 DeferredLogger& deferred_logger)
1460 {
1461 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1462
1463 BVectorWell xw(1);
1464 xw[0].resize(this->primary_variables_.numWellEq());
1465
1466 this->linSys_.recoverSolutionWell(x, xw);
1467 updateWellState(simulator, xw, well_state, deferred_logger);
1468 }
1469
1470
1471
1472
1473 template<typename TypeTag>
1474 void
1476 computeWellRatesWithBhp(const Simulator& simulator,
1477 const Scalar& bhp,
1478 std::vector<Scalar>& well_flux,
1479 DeferredLogger& deferred_logger) const
1480 {
1481 OPM_TIMEFUNCTION();
1482 const int np = this->number_of_phases_;
1483 well_flux.resize(np, 0.0);
1484
1485 const bool allow_cf = this->getAllowCrossFlow();
1486
1487 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
1488 const int cell_idx = this->well_cells_[perf];
1489 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1490 // flux for each perforation
1491 std::vector<Scalar> mob(this->num_components_, 0.);
1492 getMobility(simulator, perf, mob, deferred_logger);
1493 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
1494 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1495 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
1496
1497 std::vector<Scalar> cq_s(this->num_components_, 0.);
1498 PerforationRates<Scalar> perf_rates;
1499 computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
1500 cq_s, perf_rates, deferred_logger);
1501
1502 for(int p = 0; p < np; ++p) {
1503 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
1504 }
1505
1506 // the solvent contribution is added to the gas potentials
1507 if constexpr (has_solvent) {
1508 const auto& pu = this->phaseUsage();
1509 assert(pu.phase_used[Gas]);
1510 const int gas_pos = pu.phase_pos[Gas];
1511 well_flux[gas_pos] += cq_s[Indices::contiSolventEqIdx];
1512 }
1513 }
1514 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1515 }
1516
1517
1518
1519 template<typename TypeTag>
1520 void
1521 StandardWell<TypeTag>::
1522 computeWellRatesWithBhpIterations(const Simulator& simulator,
1523 const Scalar& bhp,
1524 std::vector<Scalar>& well_flux,
1525 DeferredLogger& deferred_logger) const
1526 {
1527 // creating a copy of the well itself, to avoid messing up the explicit information
1528 // during this copy, the only information not copied properly is the well controls
1529 StandardWell<TypeTag> well_copy(*this);
1530 well_copy.resetDampening();
1531
1532 // iterate to get a more accurate well density
1533 // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
1534 // to replace the original one
1535 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
1536 const auto& group_state = simulator.problem().wellModel().groupState();
1537
1538 // Get the current controls.
1539 const auto& summary_state = simulator.vanguard().summaryState();
1540 auto inj_controls = well_copy.well_ecl_.isInjector()
1541 ? well_copy.well_ecl_.injectionControls(summary_state)
1542 : Well::InjectionControls(0);
1543 auto prod_controls = well_copy.well_ecl_.isProducer()
1544 ? well_copy.well_ecl_.productionControls(summary_state) :
1545 Well::ProductionControls(0);
1546
1547 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
1548 auto& ws = well_state_copy.well(this->index_of_well_);
1549 if (well_copy.well_ecl_.isInjector()) {
1550 inj_controls.bhp_limit = bhp;
1551 ws.injection_cmode = Well::InjectorCMode::BHP;
1552 } else {
1553 prod_controls.bhp_limit = bhp;
1554 ws.production_cmode = Well::ProducerCMode::BHP;
1555 }
1556 ws.bhp = bhp;
1557
1558 // initialized the well rates with the potentials i.e. the well rates based on bhp
1559 const int np = this->number_of_phases_;
1560 const Scalar sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1561 for (int phase = 0; phase < np; ++phase){
1562 well_state_copy.wellRates(this->index_of_well_)[phase]
1563 = sign * ws.well_potentials[phase];
1564 }
1565 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1566 well_copy.computeAccumWell();
1567
1568 const double dt = simulator.timeStepSize();
1569 const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1570 if (!converged) {
1571 const std::string msg = " well " + name() + " did not get converged during well potential calculations "
1572 " potentials are computed based on unconverged solution";
1573 deferred_logger.debug(msg);
1574 }
1575 well_copy.updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
1576 well_copy.computeWellConnectionPressures(simulator, well_state_copy, deferred_logger);
1577 well_copy.computeWellRatesWithBhp(simulator, bhp, well_flux, deferred_logger);
1578 }
1579
1580
1581
1582
1583 template<typename TypeTag>
1584 std::vector<typename StandardWell<TypeTag>::Scalar>
1585 StandardWell<TypeTag>::
1586 computeWellPotentialWithTHP(const Simulator& simulator,
1587 DeferredLogger& deferred_logger,
1588 const WellState<Scalar>& well_state) const
1589 {
1590 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
1591 const auto& summary_state = simulator.vanguard().summaryState();
1592
1593 const auto& well = this->well_ecl_;
1594 if (well.isInjector()){
1595 const auto& controls = this->well_ecl_.injectionControls(summary_state);
1596 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
1597 if (bhp_at_thp_limit) {
1598 const Scalar bhp = std::min(*bhp_at_thp_limit,
1599 static_cast<Scalar>(controls.bhp_limit));
1600 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1601 } else {
1602 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1603 "Failed in getting converged thp based potential calculation for well "
1604 + name() + ". Instead the bhp based value is used");
1605 const Scalar bhp = controls.bhp_limit;
1606 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1607 }
1608 } else {
1609 computeWellRatesWithThpAlqProd(
1610 simulator, summary_state,
1611 deferred_logger, potentials, this->getALQ(well_state)
1612 );
1613 }
1614
1615 return potentials;
1616 }
1617
1618 template<typename TypeTag>
1619 bool
1620 StandardWell<TypeTag>::
1621 computeWellPotentialsImplicit(const Simulator& simulator,
1622 const WellState<Scalar>& well_state,
1623 std::vector<Scalar>& well_potentials,
1624 DeferredLogger& deferred_logger) const
1625 {
1626 // Create a copy of the well.
1627 // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials
1628 // is allready a copy, but not from other calls.
1629 StandardWell<TypeTag> well_copy(*this);
1630
1631 // store a copy of the well state, we don't want to update the real well state
1632 WellState<Scalar> well_state_copy = well_state;
1633 const auto& group_state = simulator.problem().wellModel().groupState();
1634 auto& ws = well_state_copy.well(this->index_of_well_);
1635
1636 // get current controls
1637 const auto& summary_state = simulator.vanguard().summaryState();
1638 auto inj_controls = well_copy.well_ecl_.isInjector()
1639 ? well_copy.well_ecl_.injectionControls(summary_state)
1640 : Well::InjectionControls(0);
1641 auto prod_controls = well_copy.well_ecl_.isProducer()
1642 ? well_copy.well_ecl_.productionControls(summary_state) :
1643 Well::ProductionControls(0);
1644
1645 // prepare/modify well state and control
1646 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
1647
1648 // update connection pressures relative to updated bhp to get better estimate of connection dp
1649 const int num_perf = ws.perf_data.size();
1650 for (int perf = 0; perf < num_perf; ++perf) {
1651 ws.perf_data.pressure[perf] = ws.bhp + well_copy.connections_.pressure_diff(perf);
1652 }
1653 // initialize rates from previous potentials
1654 const int np = this->number_of_phases_;
1655 bool trivial = true;
1656 for (int phase = 0; phase < np; ++phase){
1657 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
1658 }
1659 if (!trivial) {
1660 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
1661 for (int phase = 0; phase < np; ++phase) {
1662 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
1663 }
1664 }
1665
1666 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
1667 const double dt = simulator.timeStepSize();
1668 // iterate to get a solution at the given bhp.
1669 bool converged = false;
1670 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
1671 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1672 } else {
1673 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
1674 }
1675
1676 // fetch potentials (sign is updated on the outside).
1677 well_potentials.clear();
1678 well_potentials.resize(np, 0.0);
1679 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1680 if (has_solvent && comp_idx == Indices::contiSolventEqIdx) continue; // we do not store the solvent in the well_potentials
1681 const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
1682 well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
1683 }
1684
1685 // the solvent contribution is added to the gas potentials
1686 if constexpr (has_solvent) {
1687 const auto& pu = this->phaseUsage();
1688 assert(pu.phase_used[Gas]);
1689 const int gas_pos = pu.phase_pos[Gas];
1690 const EvalWell rate = well_copy.primary_variables_.getQs(Indices::contiSolventEqIdx);
1691 well_potentials[gas_pos] += rate.value();
1692 }
1693 return converged;
1694 }
1695
1696
1697 template<typename TypeTag>
1698 typename StandardWell<TypeTag>::Scalar
1699 StandardWell<TypeTag>::
1700 computeWellRatesAndBhpWithThpAlqProd(const Simulator &simulator,
1701 const SummaryState &summary_state,
1702 DeferredLogger& deferred_logger,
1703 std::vector<Scalar>& potentials,
1704 Scalar alq) const
1705 {
1706 Scalar bhp;
1707 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1708 simulator, summary_state, alq, deferred_logger, /*iterate_if_no_solution */ true);
1709 if (bhp_at_thp_limit) {
1710 const auto& controls = this->well_ecl_.productionControls(summary_state);
1711 bhp = std::max(*bhp_at_thp_limit,
1712 static_cast<Scalar>(controls.bhp_limit));
1713 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1714 }
1715 else {
1716 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1717 "Failed in getting converged thp based potential calculation for well "
1718 + name() + ". Instead the bhp based value is used");
1719 const auto& controls = this->well_ecl_.productionControls(summary_state);
1720 bhp = controls.bhp_limit;
1721 computeWellRatesWithBhp(simulator, bhp, potentials, deferred_logger);
1722 }
1723 return bhp;
1724 }
1725
1726 template<typename TypeTag>
1727 void
1728 StandardWell<TypeTag>::
1729 computeWellRatesWithThpAlqProd(const Simulator& simulator,
1730 const SummaryState& summary_state,
1731 DeferredLogger& deferred_logger,
1732 std::vector<Scalar>& potentials,
1733 Scalar alq) const
1734 {
1735 /*double bhp =*/
1736 computeWellRatesAndBhpWithThpAlqProd(simulator,
1737 summary_state,
1738 deferred_logger,
1739 potentials,
1740 alq);
1741 }
1742
1743 template<typename TypeTag>
1744 void
1746 computeWellPotentials(const Simulator& simulator,
1747 const WellState<Scalar>& well_state,
1748 std::vector<Scalar>& well_potentials,
1749 DeferredLogger& deferred_logger) // const
1750 {
1751 const auto [compute_potential, bhp_controlled_well] =
1752 this->WellInterfaceGeneric<Scalar>::computeWellPotentials(well_potentials, well_state);
1753
1754 if (!compute_potential) {
1755 return;
1756 }
1757
1758 bool converged_implicit = false;
1759 // for newly opened wells we dont compute the potentials implicit
1760 // group controlled wells with defaulted guiderates will have zero targets as
1761 // the potentials are used to compute the well fractions.
1762 if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
1763 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
1764 }
1765 if (!converged_implicit) {
1766 // does the well have a THP related constraint?
1767 const auto& summaryState = simulator.vanguard().summaryState();
1768 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
1769 // get the bhp value based on the bhp constraints
1770 Scalar bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1771
1772 // In some very special cases the bhp pressure target are
1773 // temporary violated. This may lead to too small or negative potentials
1774 // that could lead to premature shutting of wells.
1775 // As a remedy the bhp that gives the largest potential is used.
1776 // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp,
1777 // and the potentials will be computed using the limit as expected.
1778 const auto& ws = well_state.well(this->index_of_well_);
1779 if (this->isInjector())
1780 bhp = std::max(ws.bhp, bhp);
1781 else
1782 bhp = std::min(ws.bhp, bhp);
1783
1784 assert(std::abs(bhp) != std::numeric_limits<Scalar>::max());
1785 computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger);
1786 } else {
1787 // the well has a THP related constraint
1788 well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state);
1789 }
1790 }
1791
1792 this->checkNegativeWellPotentials(well_potentials,
1793 this->param_.check_well_operability_,
1794 deferred_logger);
1795 }
1796
1797
1798
1799
1800
1801
1802
1803 template<typename TypeTag>
1804 typename StandardWell<TypeTag>::Scalar
1806 connectionDensity([[maybe_unused]] const int globalConnIdx,
1807 const int openConnIdx) const
1808 {
1809 return (openConnIdx < 0)
1810 ? 0.0
1811 : this->connections_.rho(openConnIdx);
1812 }
1813
1814
1815
1816
1817
1818 template<typename TypeTag>
1819 void
1820 StandardWell<TypeTag>::
1821 updatePrimaryVariables(const Simulator& simulator,
1822 const WellState<Scalar>& well_state,
1823 DeferredLogger& deferred_logger)
1824 {
1825 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1826
1827 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1828 this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
1829
1830 // other primary variables related to polymer injection
1831 if constexpr (Base::has_polymermw) {
1832 this->primary_variables_.updatePolyMW(well_state);
1833 }
1834
1835 this->primary_variables_.checkFinite(deferred_logger);
1836 }
1837
1838
1839
1840
1841 template<typename TypeTag>
1842 typename StandardWell<TypeTag>::Scalar
1843 StandardWell<TypeTag>::
1844 getRefDensity() const
1845 {
1846 return this->connections_.rho();
1847 }
1848
1849
1850
1851
1852 template<typename TypeTag>
1853 void
1854 StandardWell<TypeTag>::
1855 updateWaterMobilityWithPolymer(const Simulator& simulator,
1856 const int perf,
1857 std::vector<EvalWell>& mob,
1858 DeferredLogger& deferred_logger) const
1859 {
1860 const int cell_idx = this->well_cells_[perf];
1861 const auto& int_quant = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1862 const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
1863
1864 // TODO: not sure should based on the well type or injecting/producing peforations
1865 // it can be different for crossflow
1866 if (this->isInjector()) {
1867 // assume fully mixing within injecting wellbore
1868 const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
1869 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1870 mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) );
1871 }
1872
1873 if (PolymerModule::hasPlyshlog()) {
1874 // we do not calculate the shear effects for injection wells when they do not
1875 // inject polymer.
1876 if (this->isInjector() && this->wpolymer() == 0.) {
1877 return;
1878 }
1879 // compute the well water velocity with out shear effects.
1880 // TODO: do we need to turn on crossflow here?
1881 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1882 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
1883
1884 std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
1885 PerforationRates<Scalar> perf_rates;
1886 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
1887 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1888 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
1889 computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
1890 perf_rates, deferred_logger);
1891 // TODO: make area a member
1892 const Scalar area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
1893 const auto& material_law_manager = simulator.problem().materialLawManager();
1894 const auto& scaled_drainage_info =
1895 material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
1896 const Scalar swcr = scaled_drainage_info.Swcr;
1897 const EvalWell poro = this->extendEval(int_quant.porosity());
1898 const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
1899 // guard against zero porosity and no water
1900 const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
1901 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1902 EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
1903
1904 if (PolymerModule::hasShrate()) {
1905 // the equation for the water velocity conversion for the wells and reservoir are from different version
1906 // of implementation. It can be changed to be more consistent when possible.
1907 water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
1908 }
1909 const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
1910 int_quant.pvtRegionIndex(),
1911 water_velocity);
1912 // modify the mobility with the shear factor.
1913 mob[waterCompIdx] /= shear_factor;
1914 }
1915 }
1916
1917 template<typename TypeTag>
1918 void
1919 StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian) const
1920 {
1921 this->linSys_.extract(jacobian);
1922 }
1923
1924
1925 template <typename TypeTag>
1926 void
1927 StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
1928 const BVector& weights,
1929 const int pressureVarIndex,
1930 const bool use_well_weights,
1931 const WellState<Scalar>& well_state) const
1932 {
1933 this->linSys_.extractCPRPressureMatrix(jacobian,
1934 weights,
1935 pressureVarIndex,
1936 use_well_weights,
1937 *this,
1938 Bhp,
1939 well_state);
1940 }
1941
1942
1943
1944 template<typename TypeTag>
1945 typename StandardWell<TypeTag>::EvalWell
1946 StandardWell<TypeTag>::
1947 pskinwater(const Scalar throughput,
1948 const EvalWell& water_velocity,
1949 DeferredLogger& deferred_logger) const
1950 {
1951 if constexpr (Base::has_polymermw) {
1952 const int water_table_id = this->polymerWaterTable_();
1953 if (water_table_id <= 0) {
1954 OPM_DEFLOG_THROW(std::runtime_error,
1955 fmt::format("Unused SKPRWAT table id used for well {}", name()),
1956 deferred_logger);
1957 }
1958 const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
1959 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1960 // the skin pressure when injecting water, which also means the polymer concentration is zero
1961 EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
1962 pskin_water = water_table_func.eval(throughput_eval, water_velocity);
1963 return pskin_water;
1964 } else {
1965 OPM_DEFLOG_THROW(std::runtime_error,
1966 fmt::format("Polymermw is not activated, while injecting "
1967 "skin pressure is requested for well {}", name()),
1968 deferred_logger);
1969 }
1970 }
1971
1972
1973
1974
1975
1976 template<typename TypeTag>
1977 typename StandardWell<TypeTag>::EvalWell
1978 StandardWell<TypeTag>::
1979 pskin(const Scalar throughput,
1980 const EvalWell& water_velocity,
1981 const EvalWell& poly_inj_conc,
1982 DeferredLogger& deferred_logger) const
1983 {
1984 if constexpr (Base::has_polymermw) {
1985 const Scalar sign = water_velocity >= 0. ? 1.0 : -1.0;
1986 const EvalWell water_velocity_abs = abs(water_velocity);
1987 if (poly_inj_conc == 0.) {
1988 return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
1989 }
1990 const int polymer_table_id = this->polymerTable_();
1991 if (polymer_table_id <= 0) {
1992 OPM_DEFLOG_THROW(std::runtime_error,
1993 fmt::format("Unavailable SKPRPOLY table id used for well {}", name()),
1994 deferred_logger);
1995 }
1996 const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
1997 const Scalar reference_concentration = skprpolytable.refConcentration;
1998 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
1999 // the skin pressure when injecting water, which also means the polymer concentration is zero
2000 EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2001 pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2002 if (poly_inj_conc == reference_concentration) {
2003 return sign * pskin_poly;
2004 }
2005 // poly_inj_conc != reference concentration of the table, then some interpolation will be required
2006 const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2007 const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2008 return sign * pskin;
2009 } else {
2010 OPM_DEFLOG_THROW(std::runtime_error,
2011 fmt::format("Polymermw is not activated, while injecting "
2012 "skin pressure is requested for well {}", name()),
2013 deferred_logger);
2014 }
2015 }
2016
2017
2018
2019
2020
2021 template<typename TypeTag>
2022 typename StandardWell<TypeTag>::EvalWell
2023 StandardWell<TypeTag>::
2024 wpolymermw(const Scalar throughput,
2025 const EvalWell& water_velocity,
2026 DeferredLogger& deferred_logger) const
2027 {
2028 if constexpr (Base::has_polymermw) {
2029 const int table_id = this->polymerInjTable_();
2030 const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2031 const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
2032 EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2033 if (this->wpolymer() == 0.) { // not injecting polymer
2034 return molecular_weight;
2035 }
2036 molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2037 return molecular_weight;
2038 } else {
2039 OPM_DEFLOG_THROW(std::runtime_error,
2040 fmt::format("Polymermw is not activated, while injecting "
2041 "polymer molecular weight is requested for well {}", name()),
2042 deferred_logger);
2043 }
2044 }
2045
2046
2047
2048
2049
2050 template<typename TypeTag>
2051 void
2052 StandardWell<TypeTag>::
2053 updateWaterThroughput([[maybe_unused]] const double dt,
2054 WellState<Scalar>& well_state) const
2055 {
2056 if constexpr (Base::has_polymermw) {
2057 if (!this->isInjector()) {
2058 return;
2059 }
2060
2061 auto& perf_water_throughput = well_state.well(this->index_of_well_)
2062 .perf_data.water_throughput;
2063
2064 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2065 const Scalar perf_water_vel =
2066 this->primary_variables_.value(Bhp + 1 + perf);
2067
2068 // we do not consider the formation damage due to water
2069 // flowing from reservoir into wellbore
2070 if (perf_water_vel > Scalar{0}) {
2071 perf_water_throughput[perf] += perf_water_vel * dt;
2072 }
2073 }
2074 }
2075 }
2076
2077
2078
2079
2080
2081 template<typename TypeTag>
2082 void
2083 StandardWell<TypeTag>::
2084 handleInjectivityRate(const Simulator& simulator,
2085 const int perf,
2086 std::vector<EvalWell>& cq_s) const
2087 {
2088 const int cell_idx = this->well_cells_[perf];
2089 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2090 const auto& fs = int_quants.fluidState();
2091 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2092 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2093 const int wat_vel_index = Bhp + 1 + perf;
2094 const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2095
2096 // water rate is update to use the form from water velocity, since water velocity is
2097 // a primary variable now
2098 cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
2099 }
2100
2101
2102
2103
2104 template<typename TypeTag>
2105 void
2106 StandardWell<TypeTag>::
2107 handleInjectivityEquations(const Simulator& simulator,
2108 const WellState<Scalar>& well_state,
2109 const int perf,
2110 const EvalWell& water_flux_s,
2111 DeferredLogger& deferred_logger)
2112 {
2113 const int cell_idx = this->well_cells_[perf];
2114 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2115 const auto& fs = int_quants.fluidState();
2116 const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2117 const EvalWell water_flux_r = water_flux_s / b_w;
2118 const Scalar area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2119 const EvalWell water_velocity = water_flux_r / area;
2120 const int wat_vel_index = Bhp + 1 + perf;
2121
2122 // equation for the water velocity
2123 const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
2124
2125 const auto& ws = well_state.well(this->index_of_well_);
2126 const auto& perf_data = ws.perf_data;
2127 const auto& perf_water_throughput = perf_data.water_throughput;
2128 const Scalar throughput = perf_water_throughput[perf];
2129 const int pskin_index = Bhp + 1 + this->number_of_local_perforations_ + perf;
2130
2131 EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
2132 poly_conc.setValue(this->wpolymer());
2133
2134 // equation for the skin pressure
2135 const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
2136 - pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
2137
2138 StandardWellAssemble<FluidSystem,Indices>(*this).
2139 assembleInjectivityEq(eq_pskin,
2140 eq_wat_vel,
2141 pskin_index,
2142 wat_vel_index,
2143 perf,
2144 this->primary_variables_.numWellEq(),
2145 this->linSys_);
2146 }
2147
2148
2149
2150
2151
2152 template<typename TypeTag>
2153 void
2154 StandardWell<TypeTag>::
2155 checkConvergenceExtraEqs(const std::vector<Scalar>& res,
2156 ConvergenceReport& report) const
2157 {
2158 // if different types of extra equations are involved, this function needs to be refactored further
2159
2160 // checking the convergence of the extra equations related to polymer injectivity
2161 if constexpr (Base::has_polymermw) {
2162 WellConvergence(*this).
2163 checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
2164 }
2165 }
2166
2167
2168
2169
2170
2171 template<typename TypeTag>
2172 void
2173 StandardWell<TypeTag>::
2174 updateConnectionRatePolyMW(const EvalWell& cq_s_poly,
2175 const IntensiveQuantities& int_quants,
2176 const WellState<Scalar>& well_state,
2177 const int perf,
2178 std::vector<RateVector>& connectionRates,
2179 DeferredLogger& deferred_logger) const
2180 {
2181 // the source term related to transport of molecular weight
2182 EvalWell cq_s_polymw = cq_s_poly;
2183 if (this->isInjector()) {
2184 const int wat_vel_index = Bhp + 1 + perf;
2185 const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
2186 if (water_velocity > 0.) { // injecting
2187 const auto& ws = well_state.well(this->index_of_well_);
2188 const auto& perf_water_throughput = ws.perf_data.water_throughput;
2189 const Scalar throughput = perf_water_throughput[perf];
2190 const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2191 cq_s_polymw *= molecular_weight;
2192 } else {
2193 // we do not consider the molecular weight from the polymer
2194 // going-back to the wellbore through injector
2195 cq_s_polymw *= 0.;
2196 }
2197 } else if (this->isProducer()) {
2198 if (cq_s_polymw < 0.) {
2199 cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2200 } else {
2201 // we do not consider the molecular weight from the polymer
2202 // re-injecting back through producer
2203 cq_s_polymw *= 0.;
2204 }
2205 }
2206 connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2207 }
2208
2209
2210
2211
2212
2213
2214 template<typename TypeTag>
2215 std::optional<typename StandardWell<TypeTag>::Scalar>
2216 StandardWell<TypeTag>::
2217 computeBhpAtThpLimitProd(const WellState<Scalar>& well_state,
2218 const Simulator& simulator,
2219 const SummaryState& summary_state,
2220 DeferredLogger& deferred_logger) const
2221 {
2222 return computeBhpAtThpLimitProdWithAlq(simulator,
2223 summary_state,
2224 this->getALQ(well_state),
2225 deferred_logger,
2226 /*iterate_if_no_solution */ true);
2227 }
2228
2229 template<typename TypeTag>
2230 std::optional<typename StandardWell<TypeTag>::Scalar>
2231 StandardWell<TypeTag>::
2232 computeBhpAtThpLimitProdWithAlq(const Simulator& simulator,
2233 const SummaryState& summary_state,
2234 const Scalar alq_value,
2235 DeferredLogger& deferred_logger,
2236 bool iterate_if_no_solution) const
2237 {
2238 OPM_TIMEFUNCTION();
2239 // Make the frates() function.
2240 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2241 // Not solving the well equations here, which means we are
2242 // calculating at the current Fg/Fw values of the
2243 // well. This does not matter unless the well is
2244 // crossflowing, and then it is likely still a good
2245 // approximation.
2246 std::vector<Scalar> rates(3);
2247 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2248 this->adaptRatesForVFP(rates);
2249 return rates;
2250 };
2251
2252 Scalar max_pressure = 0.0;
2253 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2254 const int cell_idx = this->well_cells_[perf];
2255 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2256 const auto& fs = int_quants.fluidState();
2257 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2258 max_pressure = std::max(max_pressure, pressure_cell);
2259 }
2260 auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
2261 summary_state,
2262 max_pressure,
2263 this->connections_.rho(),
2264 alq_value,
2265 this->getTHPConstraint(summary_state),
2266 deferred_logger);
2267
2268 if (bhpAtLimit) {
2269 auto v = frates(*bhpAtLimit);
2270 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2271 return bhpAtLimit;
2272 }
2273 }
2274
2275 if (!iterate_if_no_solution)
2276 return std::nullopt;
2277
2278 auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
2279 // Solver the well iterations to see if we are
2280 // able to get a solution with an update
2281 // solution
2282 std::vector<Scalar> rates(3);
2283 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2284 this->adaptRatesForVFP(rates);
2285 return rates;
2286 };
2287
2288 bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
2289 summary_state,
2290 max_pressure,
2291 this->connections_.rho(),
2292 alq_value,
2293 this->getTHPConstraint(summary_state),
2294 deferred_logger);
2295
2296
2297 if (bhpAtLimit) {
2298 // should we use fratesIter here since fratesIter is used in computeBhpAtThpLimitProd above?
2299 auto v = frates(*bhpAtLimit);
2300 if (std::all_of(v.cbegin(), v.cend(), [](Scalar i){ return i <= 0; }) ) {
2301 return bhpAtLimit;
2302 }
2303 }
2304
2305 // we still don't get a valied solution.
2306 return std::nullopt;
2307 }
2308
2309
2310
2311 template<typename TypeTag>
2312 std::optional<typename StandardWell<TypeTag>::Scalar>
2313 StandardWell<TypeTag>::
2314 computeBhpAtThpLimitInj(const Simulator& simulator,
2315 const SummaryState& summary_state,
2316 DeferredLogger& deferred_logger) const
2317 {
2318 // Make the frates() function.
2319 auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
2320 // Not solving the well equations here, which means we are
2321 // calculating at the current Fg/Fw values of the
2322 // well. This does not matter unless the well is
2323 // crossflowing, and then it is likely still a good
2324 // approximation.
2325 std::vector<Scalar> rates(3);
2326 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2327 return rates;
2328 };
2329
2330 return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
2331 summary_state,
2332 this->connections_.rho(),
2333 1e-6,
2334 50,
2335 true,
2336 deferred_logger);
2337 }
2338
2339
2340
2341
2342
2343 template<typename TypeTag>
2344 bool
2345 StandardWell<TypeTag>::
2346 iterateWellEqWithControl(const Simulator& simulator,
2347 const double dt,
2348 const Well::InjectionControls& inj_controls,
2349 const Well::ProductionControls& prod_controls,
2350 WellState<Scalar>& well_state,
2351 const GroupState<Scalar>& group_state,
2352 DeferredLogger& deferred_logger)
2353 {
2354 const int max_iter = this->param_.max_inner_iter_wells_;
2355 int it = 0;
2356 bool converged;
2357 bool relax_convergence = false;
2358 this->regularize_ = false;
2359 do {
2360 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2361
2362 if (it > this->param_.strict_inner_iter_wells_) {
2363 relax_convergence = true;
2364 this->regularize_ = true;
2365 }
2366
2367 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2368
2369 converged = report.converged();
2370 if (converged) {
2371 break;
2372 }
2373
2374 ++it;
2375 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2376
2377 // TODO: when this function is used for well testing purposes, will need to check the controls, so that we will obtain convergence
2378 // under the most restrictive control. Based on this converged results, we can check whether to re-open the well. Either we refactor
2379 // this function or we use different functions for the well testing purposes.
2380 // We don't allow for switching well controls while computing well potentials and testing wells
2381 // updateWellControl(simulator, well_state, deferred_logger);
2382 } while (it < max_iter);
2383
2384 return converged;
2385 }
2386
2387
2388 template<typename TypeTag>
2389 bool
2390 StandardWell<TypeTag>::
2391 iterateWellEqWithSwitching(const Simulator& simulator,
2392 const double dt,
2393 const Well::InjectionControls& inj_controls,
2394 const Well::ProductionControls& prod_controls,
2395 WellState<Scalar>& well_state,
2396 const GroupState<Scalar>& group_state,
2397 DeferredLogger& deferred_logger,
2398 const bool fixed_control /*false*/,
2399 const bool fixed_status /*false*/)
2400 {
2401 const int max_iter = this->param_.max_inner_iter_wells_;
2402 int it = 0;
2403 bool converged = false;
2404 bool relax_convergence = false;
2405 this->regularize_ = false;
2406 const auto& summary_state = simulator.vanguard().summaryState();
2407
2408 // Always take a few (more than one) iterations after a switch before allowing a new switch
2409 // The optimal number here is subject to further investigation, but it has been observerved
2410 // that unless this number is >1, we may get stuck in a cycle
2411 constexpr int min_its_after_switch = 4;
2412 int its_since_last_switch = min_its_after_switch;
2413 int switch_count= 0;
2414 // if we fail to solve eqs, we reset status/operability before leaving
2415 const auto well_status_orig = this->wellStatus_;
2416 const auto operability_orig = this->operability_status_;
2417 auto well_status_cur = well_status_orig;
2418 int status_switch_count = 0;
2419 // don't allow opening wells that are stopped from schedule or has a stopped well state
2420 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
2421 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
2422 // don't allow switcing for wells under zero rate target or requested fixed status and control
2423 const bool allow_switching =
2424 !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
2425 (!fixed_control || !fixed_status) && allow_open;
2426
2427 bool changed = false;
2428 bool final_check = false;
2429 // well needs to be set operable or else solving/updating of re-opened wells is skipped
2430 this->operability_status_.resetOperability();
2431 this->operability_status_.solvable = true;
2432 do {
2433 its_since_last_switch++;
2434 if (allow_switching && its_since_last_switch >= min_its_after_switch){
2435 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
2436 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
2437 inj_controls, prod_controls, wqTotal,
2438 deferred_logger, fixed_control, fixed_status);
2439 if (changed){
2440 its_since_last_switch = 0;
2441 switch_count++;
2442 if (well_status_cur != this->wellStatus_) {
2443 well_status_cur = this->wellStatus_;
2444 status_switch_count++;
2445 }
2446 }
2447 if (!changed && final_check) {
2448 break;
2449 } else {
2450 final_check = false;
2451 }
2452 }
2453
2454 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2455
2456 if (it > this->param_.strict_inner_iter_wells_) {
2457 relax_convergence = true;
2458 this->regularize_ = true;
2459 }
2460
2461 auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
2462
2463 converged = report.converged();
2464 if (converged) {
2465 // if equations are sufficiently linear they might converge in less than min_its_after_switch
2466 // in this case, make sure all constraints are satisfied before returning
2467 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
2468 final_check = true;
2469 its_since_last_switch = min_its_after_switch;
2470 } else {
2471 break;
2472 }
2473 }
2474
2475 ++it;
2476 solveEqAndUpdateWellState(simulator, well_state, deferred_logger);
2477
2478 } while (it < max_iter);
2479
2480 if (converged) {
2481 if (allow_switching){
2482 // update operability if status change
2483 const bool is_stopped = this->wellIsStopped();
2484 if (this->wellHasTHPConstraints(summary_state)){
2485 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
2486 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
2487 } else {
2488 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
2489 }
2490 }
2491 } else {
2492 this->wellStatus_ = well_status_orig;
2493 this->operability_status_ = operability_orig;
2494 const std::string message = fmt::format(" Well {} did not converge in {} inner iterations ("
2495 "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count);
2496 deferred_logger.debug(message);
2497 // add operability here as well ?
2498 }
2499 return converged;
2500 }
2501
2502 template<typename TypeTag>
2503 std::vector<typename StandardWell<TypeTag>::Scalar>
2505 computeCurrentWellRates(const Simulator& simulator,
2506 DeferredLogger& deferred_logger) const
2507 {
2508 // Calculate the rates that follow from the current primary variables.
2509 std::vector<Scalar> well_q_s(this->num_components_, 0.);
2510 const EvalWell& bhp = this->primary_variables_.eval(Bhp);
2511 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2512 for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
2513 const int cell_idx = this->well_cells_[perf];
2514 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
2515 std::vector<Scalar> mob(this->num_components_, 0.);
2516 getMobility(simulator, perf, mob, deferred_logger);
2517 std::vector<Scalar> cq_s(this->num_components_, 0.);
2518 Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
2519 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2520 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
2521 PerforationRates<Scalar> perf_rates;
2522 computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2523 cq_s, perf_rates, deferred_logger);
2524 for (int comp = 0; comp < this->num_components_; ++comp) {
2525 well_q_s[comp] += cq_s[comp];
2526 }
2527 }
2528 const auto& comm = this->parallel_well_info_.communication();
2529 if (comm.size() > 1)
2530 {
2531 comm.sum(well_q_s.data(), well_q_s.size());
2532 }
2533 return well_q_s;
2534 }
2535
2536
2537
2538 template <typename TypeTag>
2539 std::vector<typename StandardWell<TypeTag>::Scalar>
2541 getPrimaryVars() const
2542 {
2543 const int num_pri_vars = this->primary_variables_.numWellEq();
2544 std::vector<Scalar> retval(num_pri_vars);
2545 for (int ii = 0; ii < num_pri_vars; ++ii) {
2546 retval[ii] = this->primary_variables_.value(ii);
2547 }
2548 return retval;
2549 }
2550
2551
2552
2553
2554
2555 template <typename TypeTag>
2556 int
2557 StandardWell<TypeTag>::
2558 setPrimaryVars(typename std::vector<Scalar>::const_iterator it)
2559 {
2560 const int num_pri_vars = this->primary_variables_.numWellEq();
2561 for (int ii = 0; ii < num_pri_vars; ++ii) {
2562 this->primary_variables_.setValue(ii, it[ii]);
2563 }
2564 return num_pri_vars;
2565 }
2566
2567
2568 template <typename TypeTag>
2569 typename StandardWell<TypeTag>::Eval
2570 StandardWell<TypeTag>::
2571 connectionRateEnergy(const Scalar maxOilSaturation,
2572 const std::vector<EvalWell>& cq_s,
2573 const IntensiveQuantities& intQuants,
2574 DeferredLogger& deferred_logger) const
2575 {
2576 auto fs = intQuants.fluidState();
2577 Eval result = 0;
2578 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2579 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2580 continue;
2581 }
2582
2583 // convert to reservoir conditions
2584 EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
2585 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2586 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
2587 if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) {
2588 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2589 } else {
2590 // remove dissolved gas and vapporized oil
2591 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
2592 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
2593 // q_os = q_or * b_o + rv * q_gr * b_g
2594 // q_gs = q_gr * g_g + rs * q_or * b_o
2595 // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
2596 // d = 1.0 - rs * rv
2597 const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
2598 if (d <= 0.0) {
2599 deferred_logger.debug(
2600 fmt::format("Problematic d value {} obtained for well {}"
2601 " during calculateSinglePerf with rs {}"
2602 ", rv {}. Continue as if no dissolution (rs = 0) and"
2603 " vaporization (rv = 0) for this connection.",
2604 d, this->name(), fs.Rs(), fs.Rv()));
2605 cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
2606 } else {
2607 if (FluidSystem::gasPhaseIdx == phaseIdx) {
2608 cq_r_thermal = (cq_s[gasCompIdx] -
2609 this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) /
2610 (d * this->extendEval(fs.invB(phaseIdx)) );
2611 } else if (FluidSystem::oilPhaseIdx == phaseIdx) {
2612 // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
2613 cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) *
2614 cq_s[gasCompIdx]) /
2615 (d * this->extendEval(fs.invB(phaseIdx)) );
2616 }
2617 }
2618 }
2619
2620 // change temperature for injecting fluids
2621 if (this->isInjector() && !this->wellIsStopped() && cq_r_thermal > 0.0){
2622 // only handles single phase injection now
2623 assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
2624 fs.setTemperature(this->well_ecl_.inj_temperature());
2625 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
2626 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
2627 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
2628 paramCache.setRegionIndex(pvtRegionIdx);
2629 paramCache.setMaxOilSat(maxOilSaturation);
2630 paramCache.updatePhase(fs, phaseIdx);
2631
2632 const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
2633 fs.setDensity(phaseIdx, rho);
2634 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
2635 fs.setEnthalpy(phaseIdx, h);
2636 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2637 result += getValue(cq_r_thermal);
2638 } else if (cq_r_thermal > 0.0) {
2639 cq_r_thermal *= getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx));
2640 result += Base::restrictEval(cq_r_thermal);
2641 } else {
2642 // compute the thermal flux
2643 cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
2644 result += Base::restrictEval(cq_r_thermal);
2645 }
2646 }
2647
2648 return result * this->well_efficiency_factor_;
2649 }
2650} // namespace Opm
2651
2652#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:97
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:312
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:299
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:1424
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