My Project
Loading...
Searching...
No Matches
MultisegmentWell_impl.hpp
1/*
2 Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3 Copyright 2017 Statoil ASA.
4
5 This file is part of the Open Porous Media project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#include <opm/common/Exceptions.hpp>
22#include <opm/common/OpmLog/OpmLog.hpp>
23
24#include <opm/input/eclipse/Schedule/MSW/Segment.hpp>
25#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
26#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
27#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
28#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
29
30#include <opm/input/eclipse/Units/Units.hpp>
31
32#include <opm/material/densead/EvaluationFormat.hpp>
33
34#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
35#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
36#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
37
38#include <algorithm>
39#include <cstddef>
40#include <string>
41
42#if HAVE_CUDA || HAVE_OPENCL
43#include <opm/simulators/linalg/bda/WellContributions.hpp>
44#endif
45
46namespace Opm
47{
48
49
50 template <typename TypeTag>
51 MultisegmentWell<TypeTag>::
52 MultisegmentWell(const Well& well,
53 const ParallelWellInfo& pw_info,
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>& perf_data)
62 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
63 , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
64 , regularize_(false)
65 , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
66 {
67 // not handling solvent or polymer for now with multisegment well
68 if constexpr (has_solvent) {
69 OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
70 }
71
72 if constexpr (has_polymer) {
73 OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
74 }
75
76 if constexpr (Base::has_energy) {
77 OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
78 }
79
80 if constexpr (Base::has_foam) {
81 OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet");
82 }
83
84 if constexpr (Base::has_brine) {
85 OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
86 }
87
88 if constexpr (Base::has_watVapor) {
89 OPM_THROW(std::runtime_error, "water evaporation is not supported by multisegment well yet");
90 }
91
92 if(this->rsRvInj() > 0) {
93 OPM_THROW(std::runtime_error,
94 "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
95 " \n See (WCONINJE item 10 / WCONHIST item 8)");
96 }
97 if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
98 OPM_THROW(std::runtime_error, "water + gas case not supported by multisegment well yet");
99 }
100
101 this->thp_update_iterations = true;
102 }
103
104
105
106
107
108 template <typename TypeTag>
109 void
110 MultisegmentWell<TypeTag>::
111 init(const PhaseUsage* phase_usage_arg,
112 const std::vector<double>& depth_arg,
113 const double gravity_arg,
114 const int num_cells,
115 const std::vector< Scalar >& B_avg,
116 const bool changed_to_open_this_step)
117 {
118 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
119
120 // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
121 // for MultisegmentWell, it is much more complicated.
122 // It can be specified directly, it can be calculated from the segment depth,
123 // it can also use the cell center, which is the same for StandardWell.
124 // For the last case, should we update the depth with the depth_arg? For the
125 // future, it can be a source of wrong result with Multisegment well.
126 // An indicator from the opm-parser should indicate what kind of depth we should use here.
127
128 // \Note: we do not update the depth here. And it looks like for now, we only have the option to use
129 // specified perforation depth
130 this->initMatrixAndVectors(num_cells);
131
132 // calculate the depth difference between the perforations and the perforated grid block
133 for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
134 const int cell_idx = this->well_cells_[perf];
135 this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
136 }
137 }
138
139
140
141
142
143 template <typename TypeTag>
144 void
145 MultisegmentWell<TypeTag>::
146 initPrimaryVariablesEvaluation()
147 {
148 this->primary_variables_.init();
149 }
150
151
152
153
154
155 template <typename TypeTag>
156 void
157 MultisegmentWell<TypeTag>::
158 updatePrimaryVariables(const SummaryState& summary_state,
159 const WellState& well_state,
160 DeferredLogger& /* deferred_logger */)
161 {
162 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
163 this->primary_variables_.update(well_state, stop_or_zero_rate_target);
164 }
165
166
167
168
169
170
171 template <typename TypeTag>
172 void
175 const GroupState& group_state,
176 WellState& well_state,
178 {
179 Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
180 // scale segment rates based on the wellRates
181 // and segment pressure based on bhp
182 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
183 this->segments_.perforations(),
184 well_state);
185 this->scaleSegmentPressuresWithBhp(well_state);
186 }
187
188
189
190
191
192 template <typename TypeTag>
195 getWellConvergence(const SummaryState& /* summary_state */,
196 const WellState& well_state,
197 const std::vector<double>& B_avg,
199 const bool relax_tolerance) const
200 {
201 return this->MSWEval::getWellConvergence(well_state,
202 B_avg,
204 this->param_.max_residual_allowed_,
205 this->param_.tolerance_wells_,
206 this->param_.relaxed_tolerance_flow_well_,
207 this->param_.tolerance_pressure_ms_wells_,
208 this->param_.relaxed_tolerance_pressure_ms_well_,
210
211 }
212
213
214
215
216
217 template <typename TypeTag>
218 void
220 apply(const BVector& x, BVector& Ax) const
221 {
222 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
223 return;
224 }
225
226 if (this->param_.matrix_add_well_contributions_) {
227 // Contributions are already in the matrix itself
228 return;
229 }
230
231 this->linSys_.apply(x, Ax);
232 }
233
234
235
236
237
238 template <typename TypeTag>
239 void
241 apply(BVector& r) const
242 {
243 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
244 return;
245 }
246
247 this->linSys_.apply(r);
248 }
249
250
251
252 template <typename TypeTag>
253 void
256 const BVector& x,
257 WellState& well_state,
259 {
260 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
261 return;
262 }
263
264 BVectorWell xw(1);
265 this->linSys_.recoverSolutionWell(x, xw);
266 updateWellState(summary_state, xw, well_state, deferred_logger);
267 }
268
269
270
271
272
273 template <typename TypeTag>
274 void
276 computeWellPotentials(const Simulator& ebosSimulator,
277 const WellState& well_state,
278 std::vector<double>& well_potentials,
280 {
282 this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
283
284 if (!compute_potential) {
285 return;
286 }
287
288 debug_cost_counter_ = 0;
289 // does the well have a THP related constraint?
290 const auto& summaryState = ebosSimulator.vanguard().summaryState();
291 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
292 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
293 } else {
294 well_potentials = computeWellPotentialWithTHP(
295 well_state, ebosSimulator, deferred_logger);
296 }
297 deferred_logger.debug("Cost in iterations of finding well potential for well "
298 + this->name() + ": " + std::to_string(debug_cost_counter_));
299
300 this->checkNegativeWellPotentials(well_potentials,
301 this->param_.check_well_operability_,
303 }
304
305
306
307
308 template<typename TypeTag>
309 void
311 computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
312 std::vector<double>& well_flux,
314 {
315 if (this->well_ecl_.isInjector()) {
316 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
317 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
318 } else {
319 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
320 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
321 }
322 }
323
324 template<typename TypeTag>
325 void
326 MultisegmentWell<TypeTag>::
327 computeWellRatesWithBhp(const Simulator& ebosSimulator,
328 const double& bhp,
329 std::vector<double>& well_flux,
330 DeferredLogger& deferred_logger) const
331 {
332
333 const int np = this->number_of_phases_;
334
335 well_flux.resize(np, 0.0);
336 const bool allow_cf = this->getAllowCrossFlow();
337 const int nseg = this->numberOfSegments();
338 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
339 const auto& ws = well_state.well(this->indexOfWell());
340 auto segments_copy = ws.segments;
341 segments_copy.scale_pressure(bhp);
342 const auto& segment_pressure = segments_copy.pressure;
343 for (int seg = 0; seg < nseg; ++seg) {
344 for (const int perf : this->segments_.perforations()[seg]) {
345 const int cell_idx = this->well_cells_[perf];
346 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
347 // flux for each perforation
348 std::vector<Scalar> mob(this->num_components_, 0.);
349 getMobility(ebosSimulator, perf, mob, deferred_logger);
350 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
351 const double Tw = this->well_index_[perf] * trans_mult;
352
353 const Scalar seg_pressure = segment_pressure[seg];
354 std::vector<Scalar> cq_s(this->num_components_, 0.);
355 Scalar perf_press = 0.0;
356 PerforationRates perf_rates;
357 computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
358 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
359
360 for(int p = 0; p < np; ++p) {
361 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
362 }
363 }
364 }
365 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
366 }
367
368
369 template<typename TypeTag>
370 void
371 MultisegmentWell<TypeTag>::
372 computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
373 const Scalar& bhp,
374 std::vector<double>& well_flux,
375 DeferredLogger& deferred_logger) const
376 {
377 // creating a copy of the well itself, to avoid messing up the explicit information
378 // during this copy, the only information not copied properly is the well controls
379 MultisegmentWell<TypeTag> well_copy(*this);
380 well_copy.debug_cost_counter_ = 0;
381
382 // store a copy of the well state, we don't want to update the real well state
383 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
384 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
385 auto& ws = well_state_copy.well(this->index_of_well_);
386
387 // Get the current controls.
388 const auto& summary_state = ebosSimulator.vanguard().summaryState();
389 auto inj_controls = well_copy.well_ecl_.isInjector()
390 ? well_copy.well_ecl_.injectionControls(summary_state)
391 : Well::InjectionControls(0);
392 auto prod_controls = well_copy.well_ecl_.isProducer()
393 ? well_copy.well_ecl_.productionControls(summary_state) :
394 Well::ProductionControls(0);
395
396 // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
397 if (well_copy.well_ecl_.isInjector()) {
398 inj_controls.bhp_limit = bhp;
399 ws.injection_cmode = Well::InjectorCMode::BHP;
400 } else {
401 prod_controls.bhp_limit = bhp;
402 ws.production_cmode = Well::ProducerCMode::BHP;
403 }
404 ws.bhp = bhp;
405 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
406
407 // initialized the well rates with the potentials i.e. the well rates based on bhp
408 const int np = this->number_of_phases_;
409 bool trivial = true;
410 for (int phase = 0; phase < np; ++phase){
411 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
412 }
413 if (!trivial) {
414 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
415 for (int phase = 0; phase < np; ++phase) {
416 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
417 }
418 }
419 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
420 this->segments_.perforations(),
421 well_state_copy);
422
423 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
424 const double dt = ebosSimulator.timeStepSize();
425 // iterate to get a solution at the given bhp.
426 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
427 deferred_logger);
428
429 // compute the potential and store in the flux vector.
430 well_flux.clear();
431 well_flux.resize(np, 0.0);
432 for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
433 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
434 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
435 }
436 debug_cost_counter_ += well_copy.debug_cost_counter_;
437 }
438
439
440
441 template<typename TypeTag>
442 std::vector<double>
443 MultisegmentWell<TypeTag>::
444 computeWellPotentialWithTHP(
445 const WellState& well_state,
446 const Simulator& ebos_simulator,
447 DeferredLogger& deferred_logger) const
448 {
449 std::vector<double> potentials(this->number_of_phases_, 0.0);
450 const auto& summary_state = ebos_simulator.vanguard().summaryState();
451
452 const auto& well = this->well_ecl_;
453 if (well.isInjector()){
454 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
455 if (bhp_at_thp_limit) {
456 const auto& controls = well.injectionControls(summary_state);
457 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
458 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
459 deferred_logger.debug("Converged thp based potential calculation for well "
460 + this->name() + ", at bhp = " + std::to_string(bhp));
461 } else {
462 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
463 "Failed in getting converged thp based potential calculation for well "
464 + this->name() + ". Instead the bhp based value is used");
465 const auto& controls = well.injectionControls(summary_state);
466 const double bhp = controls.bhp_limit;
467 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
468 }
469 } else {
470 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
471 well_state, ebos_simulator, summary_state, deferred_logger);
472 if (bhp_at_thp_limit) {
473 const auto& controls = well.productionControls(summary_state);
474 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
475 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
476 deferred_logger.debug("Converged thp based potential calculation for well "
477 + this->name() + ", at bhp = " + std::to_string(bhp));
478 } else {
479 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
480 "Failed in getting converged thp based potential calculation for well "
481 + this->name() + ". Instead the bhp based value is used");
482 const auto& controls = well.productionControls(summary_state);
483 const double bhp = controls.bhp_limit;
484 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
485 }
486 }
487
488 return potentials;
489 }
490
491
492
493 template <typename TypeTag>
494 void
495 MultisegmentWell<TypeTag>::
496 solveEqAndUpdateWellState(const SummaryState& summary_state,
497 WellState& well_state,
498 DeferredLogger& deferred_logger)
499 {
500 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
501
502 // We assemble the well equations, then we check the convergence,
503 // which is why we do not put the assembleWellEq here.
504 const BVectorWell dx_well = this->linSys_.solve();
505
506 updateWellState(summary_state, dx_well, well_state, deferred_logger);
507 }
508
509
510
511
512
513 template <typename TypeTag>
514 void
515 MultisegmentWell<TypeTag>::
516 computePerfCellPressDiffs(const Simulator& ebosSimulator)
517 {
518 for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
519
520 std::vector<double> kr(this->number_of_phases_, 0.0);
521 std::vector<double> density(this->number_of_phases_, 0.0);
522
523 const int cell_idx = this->well_cells_[perf];
524 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
525 const auto& fs = intQuants.fluidState();
526
527 double sum_kr = 0.;
528
529 const PhaseUsage& pu = this->phaseUsage();
530 if (pu.phase_used[Water]) {
531 const int water_pos = pu.phase_pos[Water];
532 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
533 sum_kr += kr[water_pos];
534 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
535 }
536
537 if (pu.phase_used[Oil]) {
538 const int oil_pos = pu.phase_pos[Oil];
539 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
540 sum_kr += kr[oil_pos];
541 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
542 }
543
544 if (pu.phase_used[Gas]) {
545 const int gas_pos = pu.phase_pos[Gas];
546 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
547 sum_kr += kr[gas_pos];
548 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
549 }
550
551 assert(sum_kr != 0.);
552
553 // calculate the average density
554 double average_density = 0.;
555 for (int p = 0; p < this->number_of_phases_; ++p) {
556 average_density += kr[p] * density[p];
557 }
558 average_density /= sum_kr;
559
560 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
561 }
562 }
563
564
565
566
567
568 template <typename TypeTag>
569 void
570 MultisegmentWell<TypeTag>::
571 computeInitialSegmentFluids(const Simulator& ebos_simulator)
572 {
573 for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
574 // TODO: trying to reduce the times for the surfaceVolumeFraction calculation
575 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
576 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
577 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
578 }
579 }
580 }
581
582
583
584
585
586 template <typename TypeTag>
587 void
588 MultisegmentWell<TypeTag>::
589 updateWellState(const SummaryState& summary_state,
590 const BVectorWell& dwells,
591 WellState& well_state,
592 DeferredLogger& deferred_logger,
593 const double relaxation_factor)
594 {
595 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
596
597 const double dFLimit = this->param_.dwell_fraction_max_;
598 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
599 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
600 this->primary_variables_.updateNewton(dwells,
601 relaxation_factor,
602 dFLimit,
603 stop_or_zero_rate_target,
604 max_pressure_change);
605
606 this->primary_variables_.copyToWellState(*this, getRefDensity(), stop_or_zero_rate_target,
607 well_state, summary_state, deferred_logger);
608
609 {
610 auto& ws = well_state.well(this->index_of_well_);
611 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
612 }
613
614 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
615 }
616
617
618
619
620
621 template <typename TypeTag>
622 void
623 MultisegmentWell<TypeTag>::
624 calculateExplicitQuantities(const Simulator& ebosSimulator,
625 const WellState& well_state,
626 DeferredLogger& deferred_logger)
627 {
628 const auto& summary_state = ebosSimulator.vanguard().summaryState();
629 updatePrimaryVariables(summary_state, well_state, deferred_logger);
630 initPrimaryVariablesEvaluation();
631 computePerfCellPressDiffs(ebosSimulator);
632 computeInitialSegmentFluids(ebosSimulator);
633 }
634
635
636
637
638
639 template<typename TypeTag>
640 void
641 MultisegmentWell<TypeTag>::
642 updateProductivityIndex(const Simulator& ebosSimulator,
643 const WellProdIndexCalculator& wellPICalc,
644 WellState& well_state,
645 DeferredLogger& deferred_logger) const
646 {
647 auto fluidState = [&ebosSimulator, this](const int perf)
648 {
649 const auto cell_idx = this->well_cells_[perf];
650 return ebosSimulator.model()
651 .intensiveQuantities(cell_idx, /*timeIdx=*/ 0).fluidState();
652 };
653
654 const int np = this->number_of_phases_;
655 auto setToZero = [np](double* x) -> void
656 {
657 std::fill_n(x, np, 0.0);
658 };
659
660 auto addVector = [np](const double* src, double* dest) -> void
661 {
662 std::transform(src, src + np, dest, dest, std::plus<>{});
663 };
664
665 auto& ws = well_state.well(this->index_of_well_);
666 auto& perf_data = ws.perf_data;
667 auto* connPI = perf_data.prod_index.data();
668 auto* wellPI = ws.productivity_index.data();
669
670 setToZero(wellPI);
671
672 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
673 auto subsetPerfID = 0;
674
675 for ( const auto& perf : *this->perf_data_){
676 auto allPerfID = perf.ecl_index;
677
678 auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
679 {
680 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
681 };
682
683 std::vector<Scalar> mob(this->num_components_, 0.0);
684 getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
685
686 const auto& fs = fluidState(subsetPerfID);
687 setToZero(connPI);
688
689 if (this->isInjector()) {
690 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
691 mob, connPI, deferred_logger);
692 }
693 else { // Production or zero flow rate
694 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
695 }
696
697 addVector(connPI, wellPI);
698
699 ++subsetPerfID;
700 connPI += np;
701 }
702
703 assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
704 "Internal logic error in processing connections for PI/II");
705 }
706
707
708
709
710
711 template<typename TypeTag>
712 double
713 MultisegmentWell<TypeTag>::
714 connectionDensity(const int globalConnIdx,
715 [[maybe_unused]] const int openConnIdx) const
716 {
717 // Simple approximation: Mixture density at reservoir connection is
718 // mixture density at connection's segment.
719
720 const auto segNum = this->wellEcl()
721 .getConnections()[globalConnIdx].segment();
722
723 const auto segIdx = this->wellEcl()
724 .getSegments().segmentNumberToIndex(segNum);
725
726 return this->segments_.density(segIdx).value();
727 }
728
729
730
731
732
733 template<typename TypeTag>
734 void
735 MultisegmentWell<TypeTag>::
736 addWellContributions(SparseMatrixAdapter& jacobian) const
737 {
738 this->linSys_.extract(jacobian);
739 }
740
741
742 template<typename TypeTag>
743 void
744 MultisegmentWell<TypeTag>::
745 addWellPressureEquations(PressureMatrix& jacobian,
746 const BVector& weights,
747 const int pressureVarIndex,
748 const bool use_well_weights,
749 const WellState& well_state) const
750 {
751 // Add the pressure contribution to the cpr system for the well
752 this->linSys_.extractCPRPressureMatrix(jacobian,
753 weights,
754 pressureVarIndex,
755 use_well_weights,
756 *this,
757 this->SPres,
758 well_state);
759 }
760
761
762 template<typename TypeTag>
763 template<class Value>
764 void
765 MultisegmentWell<TypeTag>::
766 computePerfRate(const Value& pressure_cell,
767 const Value& rs,
768 const Value& rv,
769 const std::vector<Value>& b_perfcells,
770 const std::vector<Value>& mob_perfcells,
771 const double Tw,
772 const int perf,
773 const Value& segment_pressure,
774 const Value& segment_density,
775 const bool& allow_cf,
776 const std::vector<Value>& cmix_s,
777 std::vector<Value>& cq_s,
778 Value& perf_press,
779 PerforationRates& perf_rates,
780 DeferredLogger& deferred_logger) const
781 {
782 // pressure difference between the segment and the perforation
783 const Value perf_seg_press_diff = this->gravity() * segment_density *
784 this->segments_.perforation_depth_diff(perf);
785 // pressure difference between the perforation and the grid cell
786 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
787
788 // perforation pressure is the wellbore pressure corrected to perforation depth
789 // (positive sign due to convention in segments_.perforation_depth_diff() )
790 perf_press = segment_pressure + perf_seg_press_diff;
791
792 // cell pressure corrected to perforation depth
793 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
794
795 // Pressure drawdown (also used to determine direction of flow)
796 const Value drawdown = cell_press_at_perf - perf_press;
797
798 // producing perforations
799 if (drawdown > 0.0) {
800 // Do nothing if crossflow is not allowed
801 if (!allow_cf && this->isInjector()) {
802 return;
803 }
804
805 // compute component volumetric rates at standard conditions
806 for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
807 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
808 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
809 }
810
811 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
812 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
813 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
814 const Value cq_s_oil = cq_s[oilCompIdx];
815 const Value cq_s_gas = cq_s[gasCompIdx];
816 cq_s[gasCompIdx] += rs * cq_s_oil;
817 cq_s[oilCompIdx] += rv * cq_s_gas;
818 }
819 } else { // injecting perforations
820 // Do nothing if crossflow is not allowed
821 if (!allow_cf && this->isProducer()) {
822 return;
823 }
824
825 // for injecting perforations, we use total mobility
826 Value total_mob = mob_perfcells[0];
827 for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
828 total_mob += mob_perfcells[comp_idx];
829 }
830
831 // injection perforations total volume rates
832 const Value cqt_i = - Tw * (total_mob * drawdown);
833
834 // compute volume ratio between connection and at standard conditions
835 Value volume_ratio = 0.0;
836 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
837 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
838 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
839 }
840
841 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
842 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
843 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
844
845 // Incorporate RS/RV factors if both oil and gas active
846 // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
847 // basically, for injecting perforations, the wellbore is the upstreaming side.
848 const Value d = 1.0 - rv * rs;
849
850 if (getValue(d) == 0.0) {
851 OPM_DEFLOG_THROW(NumericalProblem,
852 fmt::format("Zero d value obtained for well {} "
853 "during flux calculation with rs {} and rv {}",
854 this->name(), rs, rv),
855 deferred_logger);
856 }
857
858 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
859 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
860
861 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
862 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
863 } else { // not having gas and oil at the same time
864 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
865 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
866 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
867 }
868 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
869 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
870 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
871 }
872 }
873 // injecting connections total volumerates at standard conditions
874 Value cqt_is = cqt_i / volume_ratio;
875 for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
876 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
877 }
878 } // end for injection perforations
879
880 // calculating the perforation solution gas rate and solution oil rates
881 if (this->isProducer()) {
882 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
883 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
884 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
885 // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
886 // s means standard condition, r means reservoir condition
887 // q_os = q_or * b_o + rv * q_gr * b_g
888 // q_gs = q_gr * g_g + rs * q_or * b_o
889 // d = 1.0 - rs * rv
890 // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
891 // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
892
893 const double d = 1.0 - getValue(rv) * getValue(rs);
894 // vaporized oil into gas
895 // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
896 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
897 // dissolved of gas in oil
898 // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
899 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
900 }
901 }
902 }
903
904 template <typename TypeTag>
905 template<class Value>
906 void
907 MultisegmentWell<TypeTag>::
908 computePerfRate(const IntensiveQuantities& int_quants,
909 const std::vector<Value>& mob_perfcells,
910 const double Tw,
911 const int seg,
912 const int perf,
913 const Value& segment_pressure,
914 const bool& allow_cf,
915 std::vector<Value>& cq_s,
916 Value& perf_press,
917 PerforationRates& perf_rates,
918 DeferredLogger& deferred_logger) const
919
920 {
921 auto obtain = [this](const Eval& value)
922 {
923 if constexpr (std::is_same_v<Value, Scalar>) {
924 static_cast<void>(this); // suppress clang warning
925 return getValue(value);
926 } else {
927 return this->extendEval(value);
928 }
929 };
930 auto obtainN = [](const auto& value)
931 {
932 if constexpr (std::is_same_v<Value, Scalar>) {
933 return getValue(value);
934 } else {
935 return value;
936 }
937 };
938 const auto& fs = int_quants.fluidState();
939
940 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
941 const Value rs = obtain(fs.Rs());
942 const Value rv = obtain(fs.Rv());
943
944 // not using number_of_phases_ because of solvent
945 std::vector<Value> b_perfcells(this->num_components_, 0.0);
946
947 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
948 if (!FluidSystem::phaseIsActive(phaseIdx)) {
949 continue;
950 }
951
952 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
953 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
954 }
955
956 std::vector<Value> cmix_s(this->numComponents(), 0.0);
957 for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
958 cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
959 }
960
961 this->computePerfRate(pressure_cell,
962 rs,
963 rv,
964 b_perfcells,
965 mob_perfcells,
966 Tw,
967 perf,
968 segment_pressure,
969 obtainN(this->segments_.density(seg)),
970 allow_cf,
971 cmix_s,
972 cq_s,
973 perf_press,
974 perf_rates,
975 deferred_logger);
976 }
977
978 template <typename TypeTag>
979 void
980 MultisegmentWell<TypeTag>::
981 computeSegmentFluidProperties(const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
982 {
983 // TODO: the concept of phases and components are rather confusing in this function.
984 // needs to be addressed sooner or later.
985
986 // get the temperature for later use. It is only useful when we are not handling
987 // thermal related simulation
988 // basically, it is a single value for all the segments
989
990 EvalWell temperature;
991 EvalWell saltConcentration;
992 // not sure how to handle the pvt region related to segment
993 // for the current approach, we use the pvt region of the first perforated cell
994 // although there are some text indicating using the pvt region of the lowest
995 // perforated cell
996 // TODO: later to investigate how to handle the pvt region
997 int pvt_region_index;
998 {
999 // using the first perforated cell
1000 const int cell_idx = this->well_cells_[0];
1001 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1002 const auto& fs = intQuants.fluidState();
1003 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1004 saltConcentration = this->extendEval(fs.saltConcentration());
1005 pvt_region_index = fs.pvtRegionIndex();
1006 }
1007
1008 this->segments_.computeFluidProperties(temperature,
1009 saltConcentration,
1010 this->primary_variables_,
1011 pvt_region_index,
1012 deferred_logger);
1013 }
1014
1015 template <typename TypeTag>
1016 template<class Value>
1017 void
1018 MultisegmentWell<TypeTag>::
1019 getMobility(const Simulator& ebosSimulator,
1020 const int perf,
1021 std::vector<Value>& mob,
1022 DeferredLogger& deferred_logger) const
1023 {
1024 auto obtain = [this](const Eval& value)
1025 {
1026 if constexpr (std::is_same_v<Value, Scalar>) {
1027 static_cast<void>(this); // suppress clang warning
1028 return getValue(value);
1029 } else {
1030 return this->extendEval(value);
1031 }
1032 };
1033
1034 WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
1035
1036 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
1037 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1038 const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
1039 const int seg = this->segmentNumberToIndex(con.segment());
1040 // from the reference results, it looks like MSW uses segment pressure instead of BHP here
1041 // Note: this is against the documented definition.
1042 // we can change this depending on what we want
1043 const double segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
1044 const double perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
1045 * this->segments_.perforation_depth_diff(perf);
1046 const double perf_press = segment_pres + perf_seg_press_diff;
1047 const double multiplier = this->getInjMult(perf, segment_pres, perf_press);
1048 for (std::size_t i = 0; i < mob.size(); ++i) {
1049 mob[i] *= multiplier;
1050 }
1051 }
1052 }
1053
1054
1055
1056 template<typename TypeTag>
1057 double
1058 MultisegmentWell<TypeTag>::
1059 getRefDensity() const
1060 {
1061 return this->segments_.getRefDensity();
1062 }
1063
1064 template<typename TypeTag>
1065 void
1066 MultisegmentWell<TypeTag>::
1067 checkOperabilityUnderBHPLimit(const WellState& /*well_state*/, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1068 {
1069 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1070 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1071 // Crude but works: default is one atmosphere.
1072 // TODO: a better way to detect whether the BHP is defaulted or not
1073 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1074 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1075 // if the BHP limit is not defaulted or the well does not have a THP limit
1076 // we need to check the BHP limit
1077 double total_ipr_mass_rate = 0.0;
1078 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1079 {
1080 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1081 continue;
1082 }
1083
1084 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1085 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1086
1087 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1088 total_ipr_mass_rate += ipr_rate * rho;
1089 }
1090 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1091 this->operability_status_.operable_under_only_bhp_limit = false;
1092 }
1093
1094 // checking whether running under BHP limit will violate THP limit
1095 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1096 // option 1: calculate well rates based on the BHP limit.
1097 // option 2: stick with the above IPR curve
1098 // we use IPR here
1099 std::vector<double> well_rates_bhp_limit;
1100 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1101
1102 const double thp_limit = this->getTHPConstraint(summaryState);
1103 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1104 bhp_limit,
1105 this->getRefDensity(),
1106 this->wellEcl().alq_value(),
1107 thp_limit,
1108 deferred_logger);
1109 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1110 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1111 }
1112 }
1113 } else {
1114 // defaulted BHP and there is a THP constraint
1115 // default BHP limit is about 1 atm.
1116 // when applied the hydrostatic pressure correction dp,
1117 // most likely we get a negative value (bhp + dp)to search in the VFP table,
1118 // which is not desirable.
1119 // we assume we can operate under defaulted BHP limit and will violate the THP limit
1120 // when operating under defaulted BHP limit.
1121 this->operability_status_.operable_under_only_bhp_limit = true;
1122 this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1123 }
1124 }
1125
1126
1127
1128 template<typename TypeTag>
1129 void
1130 MultisegmentWell<TypeTag>::
1131 updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const
1132 {
1133 // TODO: not handling solvent related here for now
1134
1135 // initialize all the values to be zero to begin with
1136 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1137 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1138
1139 const int nseg = this->numberOfSegments();
1140 std::vector<double> seg_dp(nseg, 0.0);
1141 for (int seg = 0; seg < nseg; ++seg) {
1142 // calculating the perforation rate for each perforation that belongs to this segment
1143 const double dp = this->getSegmentDp(seg,
1144 this->segments_.density(seg).value(),
1145 seg_dp);
1146 seg_dp[seg] = dp;
1147 for (const int perf : this->segments_.perforations()[seg]) {
1148 std::vector<Scalar> mob(this->num_components_, 0.0);
1149
1150 // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
1151 getMobility(ebos_simulator, perf, mob, deferred_logger);
1152
1153 const int cell_idx = this->well_cells_[perf];
1154 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1155 const auto& fs = int_quantities.fluidState();
1156 // pressure difference between the segment and the perforation
1157 const double perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1158 // pressure difference between the perforation and the grid cell
1159 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1160 const double pressure_cell = this->getPerfCellPressure(fs).value();
1161
1162 // calculating the b for the connection
1163 std::vector<double> b_perf(this->num_components_);
1164 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1165 if (!FluidSystem::phaseIsActive(phase)) {
1166 continue;
1167 }
1168 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1169 b_perf[comp_idx] = fs.invB(phase).value();
1170 }
1171
1172 // the pressure difference between the connection and BHP
1173 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1174 const double pressure_diff = pressure_cell - h_perf;
1175
1176 // do not take into consideration the crossflow here.
1177 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1178 deferred_logger.debug("CROSSFLOW_IPR",
1179 "cross flow found when updateIPR for well " + this->name());
1180 }
1181
1182 // the well index associated with the connection
1183 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1184
1185 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1186 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1187 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1188 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1189 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1190 ipr_b_perf[comp_idx] += tw_mob;
1191 }
1192
1193 // we need to handle the rs and rv when both oil and gas are present
1194 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1195 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1196 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1197 const double rs = (fs.Rs()).value();
1198 const double rv = (fs.Rv()).value();
1199
1200 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1201 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1202
1203 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1204 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1205
1206 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1207 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1208
1209 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1210 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1211 }
1212
1213 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1214 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1215 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1216 }
1217 }
1218 }
1219 }
1220
1221 template<typename TypeTag>
1222 void
1223 MultisegmentWell<TypeTag>::
1224 checkOperabilityUnderTHPLimit(
1225 const Simulator& ebos_simulator,
1226 const WellState& well_state,
1227 DeferredLogger& deferred_logger)
1228 {
1229 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1230 const auto obtain_bhp = this->isProducer()
1231 ? computeBhpAtThpLimitProd(
1232 well_state, ebos_simulator, summaryState, deferred_logger)
1233 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1234
1235 if (obtain_bhp) {
1236 this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1237
1238 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1239 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1240
1241 const double thp_limit = this->getTHPConstraint(summaryState);
1242 if (this->isProducer() && *obtain_bhp < thp_limit) {
1243 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1244 + " bars is SMALLER than thp limit "
1245 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1246 + " bars as a producer for well " + this->name();
1247 deferred_logger.debug(msg);
1248 }
1249 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1250 const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1251 + " bars is LARGER than thp limit "
1252 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1253 + " bars as a injector for well " + this->name();
1254 deferred_logger.debug(msg);
1255 }
1256 } else {
1257 // Shutting wells that can not find bhp value from thp
1258 // when under THP control
1259 this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1260 this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1261 if (!this->wellIsStopped()) {
1262 const double thp_limit = this->getTHPConstraint(summaryState);
1263 deferred_logger.debug(" could not find bhp value at thp limit "
1264 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1265 + " bar for well " + this->name() + ", the well might need to be closed ");
1266 }
1267 }
1268 }
1269
1270
1271
1272
1273
1274 template<typename TypeTag>
1275 bool
1276 MultisegmentWell<TypeTag>::
1277 iterateWellEqWithControl(const Simulator& ebosSimulator,
1278 const double dt,
1279 const Well::InjectionControls& inj_controls,
1280 const Well::ProductionControls& prod_controls,
1281 WellState& well_state,
1282 const GroupState& group_state,
1283 DeferredLogger& deferred_logger)
1284 {
1285 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
1286
1287 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1288
1289 {
1290 // getWellFiniteResiduals returns false for nan/inf residuals
1291 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1292 if(!isFinite)
1293 return false;
1294 }
1295
1296 std::vector<std::vector<Scalar> > residual_history;
1297 std::vector<double> measure_history;
1298 int it = 0;
1299 // relaxation factor
1300 double relaxation_factor = 1.;
1301 const double min_relaxation_factor = 0.6;
1302 bool converged = false;
1303 int stagnate_count = 0;
1304 bool relax_convergence = false;
1305 this->regularize_ = false;
1306 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1307
1308 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1309
1310 const BVectorWell dx_well = this->linSys_.solve();
1311
1312 if (it > this->param_.strict_inner_iter_wells_) {
1313 relax_convergence = true;
1314 this->regularize_ = true;
1315 }
1316
1317 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1318 const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1319 if (report.converged()) {
1320 converged = true;
1321 break;
1322 }
1323
1324 {
1325 // getFinteWellResiduals returns false for nan/inf residuals
1326 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1327 if (!isFinite)
1328 return false;
1329
1330 residual_history.push_back(residuals);
1331 measure_history.push_back(this->getResidualMeasureValue(well_state,
1332 residual_history[it],
1333 this->param_.tolerance_wells_,
1334 this->param_.tolerance_pressure_ms_wells_,
1335 deferred_logger) );
1336 }
1337
1338
1339 bool is_oscillate = false;
1340 bool is_stagnate = false;
1341
1342 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1343 // TODO: maybe we should have more sophisticated strategy to recover the relaxation factor,
1344 // for example, to recover it to be bigger
1345
1346 if (is_oscillate || is_stagnate) {
1347 // HACK!
1348 std::ostringstream sstr;
1349 if (relaxation_factor == min_relaxation_factor) {
1350 // Still stagnating, terminate iterations if 5 iterations pass.
1351 ++stagnate_count;
1352 if (stagnate_count == 6) {
1353 sstr << " well " << this->name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1354 const auto reportStag = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, true);
1355 if (reportStag.converged()) {
1356 converged = true;
1357 sstr << " well " << this->name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations";
1358 deferred_logger.debug(sstr.str());
1359 return converged;
1360 }
1361 }
1362 }
1363
1364 // a factor value to reduce the relaxation_factor
1365 const double reduction_mutliplier = 0.9;
1366 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1367
1368 // debug output
1369 if (is_stagnate) {
1370 sstr << " well " << this->name() << " observes stagnation in inner iteration " << it << "\n";
1371
1372 }
1373 if (is_oscillate) {
1374 sstr << " well " << this->name() << " observes oscillation in inner iteration " << it << "\n";
1375 }
1376 sstr << " relaxation_factor is " << relaxation_factor << " now\n";
1377
1378 this->regularize_ = true;
1379 deferred_logger.debug(sstr.str());
1380 }
1381 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1382 initPrimaryVariablesEvaluation();
1383 }
1384
1385 // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
1386 if (converged) {
1387 std::ostringstream sstr;
1388 sstr << " Well " << this->name() << " converged in " << it << " inner iterations.";
1389 if (relax_convergence)
1390 sstr << " (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ << " iterations)";
1391 deferred_logger.debug(sstr.str());
1392 } else {
1393 std::ostringstream sstr;
1394 sstr << " Well " << this->name() << " did not converge in " << it << " inner iterations.";
1395#define EXTRA_DEBUG_MSW 0
1396#if EXTRA_DEBUG_MSW
1397 sstr << "***** Outputting the residual history for well " << this->name() << " during inner iterations:";
1398 for (int i = 0; i < it; ++i) {
1399 const auto& residual = residual_history[i];
1400 sstr << " residual at " << i << "th iteration ";
1401 for (const auto& res : residual) {
1402 sstr << " " << res;
1403 }
1404 sstr << " " << measure_history[i] << " \n";
1405 }
1406#endif
1407#undef EXTRA_DEBUG_MSW
1408 deferred_logger.debug(sstr.str());
1409 }
1410
1411 return converged;
1412 }
1413
1414
1415 template<typename TypeTag>
1416 bool
1417 MultisegmentWell<TypeTag>::
1418 iterateWellEqWithSwitching(const Simulator& ebosSimulator,
1419 const double dt,
1420 const Well::InjectionControls& inj_controls,
1421 const Well::ProductionControls& prod_controls,
1422 WellState& well_state,
1423 const GroupState& group_state,
1424 DeferredLogger& deferred_logger)
1425 {
1426 //if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
1427
1428 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1429
1430 {
1431 // getWellFiniteResiduals returns false for nan/inf residuals
1432 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1433 if(!isFinite)
1434 return false;
1435 }
1436
1437 std::vector<std::vector<Scalar> > residual_history;
1438 std::vector<double> measure_history;
1439 int it = 0;
1440 // relaxation factor
1441 double relaxation_factor = 1.;
1442 const double min_relaxation_factor = 0.6;
1443 bool converged = false;
1444 [[maybe_unused]] int stagnate_count = 0;
1445 bool relax_convergence = false;
1446 this->regularize_ = false;
1447
1448 // Max status switch frequency should be 2 to avoid getting stuck in cycle
1449 const int min_its_after_switch = 2;
1450 int its_since_last_switch = min_its_after_switch;
1451 int switch_count= 0;
1452 const auto well_status = this->wellStatus_;
1453 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1454 const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
1455 bool changed = false;
1456 bool final_check = false;
1457
1458 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1459 its_since_last_switch++;
1460 if (its_since_last_switch >= min_its_after_switch){
1461 const double wqTotal = this->primary_variables_.getWQTotal().value();
1462 changed = this->updateWellControlAndStatusLocalIteration (ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger);
1463 if (changed){
1464 its_since_last_switch = 0;
1465 switch_count++;
1466 }
1467 if (!changed && final_check) {
1468 break;
1469 } else {
1470 final_check = false;
1471 }
1472 }
1473
1474 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1475
1476 const BVectorWell dx_well = this->linSys_.solve();
1477
1478 if (it > this->param_.strict_inner_iter_wells_) {
1479 relax_convergence = true;
1480 this->regularize_ = true;
1481 }
1482
1483 const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1484 converged = report.converged();
1485 if (converged) {
1486 // if equations are sufficiently linear they might converge in less than min_its_after_switch
1487 // in this case, make sure all constraints are satisfied before returning
1488 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1489 final_check = true;
1490 its_since_last_switch = min_its_after_switch;
1491 } else {
1492 break;
1493 }
1494 }
1495
1496 // getFinteWellResiduals returns false for nan/inf residuals
1497 {
1498 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1499 if (!isFinite)
1500 return false;
1501
1502 residual_history.push_back(residuals);
1503 }
1504
1505 if (!converged) {
1506 measure_history.push_back(this->getResidualMeasureValue(well_state,
1507 residual_history[it],
1508 this->param_.tolerance_wells_,
1509 this->param_.tolerance_pressure_ms_wells_,
1510 deferred_logger));
1511
1512 bool is_oscillate = false;
1513 bool is_stagnate = false;
1514
1515 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1516 // TODO: maybe we should have more sophisticated strategy to recover the relaxation factor,
1517 // for example, to recover it to be bigger
1518
1519 if (is_oscillate || is_stagnate) {
1520 // HACK!
1521 std::string message;
1522 if (relaxation_factor == min_relaxation_factor) {
1523 ++stagnate_count;
1524 if (false) { // this disables the usage of the relaxed tolerance
1525 fmt::format_to(std::back_inserter(message), " Well {} observes severe stagnation and/or oscillation."
1526 " We relax the tolerance and check for convergence. \n", this->name());
1527 const auto reportStag = getWellConvergence(summary_state, well_state, Base::B_avg_,
1528 deferred_logger, true);
1529 if (reportStag.converged()) {
1530 converged = true;
1531 fmt::format_to(std::back_inserter(message), " Well {} manages to get converged with relaxed tolerances in {} inner iterations", this->name(), it);
1532 deferred_logger.debug(message);
1533 return converged;
1534 }
1535 }
1536 }
1537
1538 // a factor value to reduce the relaxation_factor
1539 constexpr double reduction_mutliplier = 0.9;
1540 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1541
1542 // debug output
1543 if (is_stagnate) {
1544 fmt::format_to(std::back_inserter(message), " well {} observes stagnation in inner iteration {}\n", this->name(), it);
1545 }
1546 if (is_oscillate) {
1547 fmt::format_to(std::back_inserter(message), " well {} observes oscillation in inner iteration {}\n", this->name(), it);
1548 }
1549 fmt::format_to(std::back_inserter(message), " relaxation_factor is {} now\n", relaxation_factor);
1550
1551 this->regularize_ = true;
1552 deferred_logger.debug(message);
1553 }
1554 }
1555 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1556 initPrimaryVariablesEvaluation();
1557 }
1558
1559 if (converged) {
1560 if (allow_switching){
1561 // update operability if status change
1562 const bool is_stopped = this->wellIsStopped();
1563 if (this->wellHasTHPConstraints(summary_state)){
1564 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
1565 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
1566 } else {
1567 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
1568 }
1569 // We reset the well status to it's original state. Status is updated
1570 // on the outside based on operability status
1571 this->wellStatus_ = well_status;
1572 }
1573 std::string message = fmt::format(" Well {} converged in {} inner iterations ("
1574 "{} control/status switches).", this->name(), it, switch_count);
1575 if (relax_convergence) {
1576 message.append(fmt::format(" (A relaxed tolerance was used after {} iterations)",
1577 this->param_.strict_inner_iter_wells_));
1578 }
1579 deferred_logger.debug(message);
1580 } else {
1581 const std::string message = fmt::format(" Well {} did not converged in {} inner iterations ("
1582 "{} control/status switches).", this->name(), it, switch_count);
1583 deferred_logger.debug(message);
1584 }
1585
1586 return converged;
1587 }
1588
1589
1590 template<typename TypeTag>
1591 void
1592 MultisegmentWell<TypeTag>::
1593 assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
1594 const double dt,
1595 const Well::InjectionControls& inj_controls,
1596 const Well::ProductionControls& prod_controls,
1597 WellState& well_state,
1598 const GroupState& group_state,
1599 DeferredLogger& deferred_logger)
1600 {
1601 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1602
1603 // update the upwinding segments
1604 this->segments_.updateUpwindingSegments(this->primary_variables_);
1605
1606 // calculate the fluid properties needed.
1607 computeSegmentFluidProperties(ebosSimulator, deferred_logger);
1608
1609 // clear all entries
1610 this->linSys_.clear();
1611
1612 auto& ws = well_state.well(this->index_of_well_);
1613 ws.phase_mixing_rates.fill(0.0);
1614
1615 // for the black oil cases, there will be four equations,
1616 // the first three of them are the mass balance equations, the last one is the pressure equations.
1617 //
1618 // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
1619
1620 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1621
1622 const int nseg = this->numberOfSegments();
1623
1624 for (int seg = 0; seg < nseg; ++seg) {
1625 // calculating the accumulation term
1626 // TODO: without considering the efficiency factor for now
1627 {
1628 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1629
1630 // Add a regularization_factor to increase the accumulation term
1631 // This will make the system less stiff and help convergence for
1632 // difficult cases
1633 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1634 // for each component
1635 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1636 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1637 - segment_fluid_initial_[seg][comp_idx]) / dt;
1638 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1639 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1640 }
1641 }
1642 // considering the contributions due to flowing out from the segment
1643 {
1644 const int seg_upwind = this->segments_.upwinding_segment(seg);
1645 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1646 const EvalWell segment_rate =
1647 this->primary_variables_.getSegmentRateUpwinding(seg,
1648 seg_upwind,
1649 comp_idx) *
1650 this->well_efficiency_factor_;
1651 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1652 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1653 }
1654 }
1655
1656 // considering the contributions from the inlet segments
1657 {
1658 for (const int inlet : this->segments_.inlets()[seg]) {
1659 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1660 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1661 const EvalWell inlet_rate =
1662 this->primary_variables_.getSegmentRateUpwinding(inlet,
1663 inlet_upwind,
1664 comp_idx) *
1665 this->well_efficiency_factor_;
1666 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1667 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1668 }
1669 }
1670 }
1671
1672 // calculating the perforation rate for each perforation that belongs to this segment
1673 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1674 auto& perf_data = ws.perf_data;
1675 auto& perf_rates = perf_data.phase_rates;
1676 auto& perf_press_state = perf_data.pressure;
1677 for (const int perf : this->segments_.perforations()[seg]) {
1678 const int cell_idx = this->well_cells_[perf];
1679 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1680 std::vector<EvalWell> mob(this->num_components_, 0.0);
1681 getMobility(ebosSimulator, perf, mob, deferred_logger);
1682 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1683 const double Tw = this->well_index_[perf] * trans_mult;
1684 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1685 EvalWell perf_press;
1686 PerforationRates perfRates;
1687 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
1688 allow_cf, cq_s, perf_press, perfRates, deferred_logger);
1689
1690 // updating the solution gas rate and solution oil rate
1691 if (this->isProducer()) {
1692 ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas;
1693 ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil;
1694 }
1695
1696 // store the perf pressure and rates
1697 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1698 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1699 }
1700 perf_press_state[perf] = perf_press.value();
1701
1702 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1703 // the cq_s entering mass balance equations need to consider the efficiency factors.
1704 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1705
1706 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1707
1708 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1709 assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
1710 }
1711 }
1712
1713 // the fourth dequation, the pressure drop equation
1714 if (seg == 0) { // top segment, pressure equation is the control equation
1715 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1716 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1717 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1718 assembleControlEq(well_state,
1719 group_state,
1720 schedule,
1721 summaryState,
1722 inj_controls,
1723 prod_controls,
1724 getRefDensity(),
1725 this->primary_variables_,
1726 this->linSys_,
1727 deferred_logger);
1728 } else {
1729 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1730 this->assemblePressureEq(seg, unit_system, well_state, this->param_.use_average_density_ms_wells_, deferred_logger);
1731 }
1732 }
1733
1734 this->linSys_.createSolver();
1735 }
1736
1737
1738
1739
1740 template<typename TypeTag>
1741 bool
1742 MultisegmentWell<TypeTag>::
1743 openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
1744 {
1745 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1746 }
1747
1748
1749 template<typename TypeTag>
1750 bool
1751 MultisegmentWell<TypeTag>::
1752 allDrawDownWrongDirection(const Simulator& ebos_simulator) const
1753 {
1754 bool all_drawdown_wrong_direction = true;
1755 const int nseg = this->numberOfSegments();
1756
1757 for (int seg = 0; seg < nseg; ++seg) {
1758 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1759 for (const int perf : this->segments_.perforations()[seg]) {
1760
1761 const int cell_idx = this->well_cells_[perf];
1762 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1763 const auto& fs = intQuants.fluidState();
1764
1765 // pressure difference between the segment and the perforation
1766 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1767 // pressure difference between the perforation and the grid cell
1768 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1769
1770 const double pressure_cell = this->getPerfCellPressure(fs).value();
1771 const double perf_press = pressure_cell - cell_perf_press_diff;
1772 // Pressure drawdown (also used to determine direction of flow)
1773 // TODO: not 100% sure about the sign of the seg_perf_press_diff
1774 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1775
1776 // for now, if there is one perforation can produce/inject in the correct
1777 // direction, we consider this well can still produce/inject.
1778 // TODO: it can be more complicated than this to cause wrong-signed rates
1779 if ( (drawdown < 0. && this->isInjector()) ||
1780 (drawdown > 0. && this->isProducer()) ) {
1781 all_drawdown_wrong_direction = false;
1782 break;
1783 }
1784 }
1785 }
1786
1787 return all_drawdown_wrong_direction;
1788 }
1789
1790
1791
1792
1793 template<typename TypeTag>
1794 void
1795 MultisegmentWell<TypeTag>::
1796 updateWaterThroughput(const double /*dt*/, WellState& /*well_state*/) const
1797 {
1798 }
1799
1800
1801
1802
1803
1804 template<typename TypeTag>
1805 typename MultisegmentWell<TypeTag>::EvalWell
1806 MultisegmentWell<TypeTag>::
1807 getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const
1808 {
1809 EvalWell temperature;
1810 EvalWell saltConcentration;
1811 int pvt_region_index;
1812 {
1813 // using the pvt region of first perforated cell
1814 // TODO: it should be a member of the WellInterface, initialized properly
1815 const int cell_idx = this->well_cells_[0];
1816 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
1817 const auto& fs = intQuants.fluidState();
1818 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1819 saltConcentration = this->extendEval(fs.saltConcentration());
1820 pvt_region_index = fs.pvtRegionIndex();
1821 }
1822
1823 return this->segments_.getSurfaceVolume(temperature,
1824 saltConcentration,
1825 this->primary_variables_,
1826 pvt_region_index,
1827 seg_idx);
1828 }
1829
1830
1831 template<typename TypeTag>
1832 std::optional<double>
1833 MultisegmentWell<TypeTag>::
1834 computeBhpAtThpLimitProd(const WellState& well_state,
1835 const Simulator& ebos_simulator,
1836 const SummaryState& summary_state,
1837 DeferredLogger& deferred_logger) const
1838 {
1839 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
1840 ebos_simulator,
1841 summary_state,
1842 this->getALQ(well_state),
1843 deferred_logger);
1844 }
1845
1846
1847
1848 template<typename TypeTag>
1849 std::optional<double>
1850 MultisegmentWell<TypeTag>::
1851 computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
1852 const SummaryState& summary_state,
1853 const double alq_value,
1854 DeferredLogger& deferred_logger) const
1855 {
1856 // Make the frates() function.
1857 auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1858 // Not solving the well equations here, which means we are
1859 // calculating at the current Fg/Fw values of the
1860 // well. This does not matter unless the well is
1861 // crossflowing, and then it is likely still a good
1862 // approximation.
1863 std::vector<double> rates(3);
1864 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1865 return rates;
1866 };
1867
1868 auto bhpAtLimit = WellBhpThpCalculator(*this).
1869 computeBhpAtThpLimitProd(frates,
1870 summary_state,
1871 this->maxPerfPress(ebos_simulator),
1872 this->getRefDensity(),
1873 alq_value,
1874 this->getTHPConstraint(summary_state),
1875 deferred_logger);
1876
1877 if (bhpAtLimit)
1878 return bhpAtLimit;
1879
1880 auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1881 // Solver the well iterations to see if we are
1882 // able to get a solution with an update
1883 // solution
1884 std::vector<double> rates(3);
1885 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1886 return rates;
1887 };
1888
1889 return WellBhpThpCalculator(*this).
1890 computeBhpAtThpLimitProd(fratesIter,
1891 summary_state,
1892 this->maxPerfPress(ebos_simulator),
1893 this->getRefDensity(),
1894 alq_value,
1895 this->getTHPConstraint(summary_state),
1896 deferred_logger);
1897 }
1898
1899 template<typename TypeTag>
1900 std::optional<double>
1901 MultisegmentWell<TypeTag>::
1902 computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
1903 const SummaryState& summary_state,
1904 DeferredLogger& deferred_logger) const
1905 {
1906 // Make the frates() function.
1907 auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1908 // Not solving the well equations here, which means we are
1909 // calculating at the current Fg/Fw values of the
1910 // well. This does not matter unless the well is
1911 // crossflowing, and then it is likely still a good
1912 // approximation.
1913 std::vector<double> rates(3);
1914 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1915 return rates;
1916 };
1917
1918 auto bhpAtLimit = WellBhpThpCalculator(*this).
1919 computeBhpAtThpLimitInj(frates,
1920 summary_state,
1921 this->getRefDensity(),
1922 0.05,
1923 100,
1924 false,
1925 deferred_logger);
1926
1927 if (bhpAtLimit)
1928 return bhpAtLimit;
1929
1930 auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1931 // Solver the well iterations to see if we are
1932 // able to get a solution with an update
1933 // solution
1934 std::vector<double> rates(3);
1935 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1936 return rates;
1937 };
1938
1939 return WellBhpThpCalculator(*this).
1940 computeBhpAtThpLimitInj(fratesIter,
1941 summary_state,
1942 this->getRefDensity(),
1943 0.05,
1944 100,
1945 false,
1946 deferred_logger);
1947 }
1948
1949
1950
1951
1952
1953 template<typename TypeTag>
1954 double
1955 MultisegmentWell<TypeTag>::
1956 maxPerfPress(const Simulator& ebos_simulator) const
1957 {
1958 double max_pressure = 0.0;
1959 const int nseg = this->numberOfSegments();
1960 for (int seg = 0; seg < nseg; ++seg) {
1961 for (const int perf : this->segments_.perforations()[seg]) {
1962 const int cell_idx = this->well_cells_[perf];
1963 const auto& int_quants = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1964 const auto& fs = int_quants.fluidState();
1965 double pressure_cell = this->getPerfCellPressure(fs).value();
1966 max_pressure = std::max(max_pressure, pressure_cell);
1967 }
1968 }
1969 return max_pressure;
1970 }
1971
1972
1973
1974
1975
1976 template<typename TypeTag>
1977 std::vector<double>
1979 computeCurrentWellRates(const Simulator& ebosSimulator,
1981 {
1982 // Calculate the rates that follow from the current primary variables.
1983 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
1984 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1985 const int nseg = this->numberOfSegments();
1986 for (int seg = 0; seg < nseg; ++seg) {
1987 // calculating the perforation rate for each perforation that belongs to this segment
1988 const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
1989 for (const int perf : this->segments_.perforations()[seg]) {
1990 const int cell_idx = this->well_cells_[perf];
1991 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
1992 std::vector<Scalar> mob(this->num_components_, 0.0);
1993 getMobility(ebosSimulator, perf, mob, deferred_logger);
1994 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1995 const double Tw = this->well_index_[perf] * trans_mult;
1996 std::vector<Scalar> cq_s(this->num_components_, 0.0);
1997 Scalar perf_press = 0.0;
1999 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
2001 for (int comp = 0; comp < this->num_components_; ++comp) {
2002 well_q_s[comp] += cq_s[comp];
2003 }
2004 }
2005 }
2006 return well_q_s;
2007 }
2008
2009
2010 template <typename TypeTag>
2011 std::vector<double>
2013 getPrimaryVars() const
2014 {
2015 const int num_seg = this->numberOfSegments();
2016 constexpr int num_eq = MSWEval::numWellEq;
2017 std::vector<double> retval(num_seg * num_eq);
2018 for (int ii = 0; ii < num_seg; ++ii) {
2019 const auto& pv = this->primary_variables_.value(ii);
2020 std::copy(pv.begin(), pv.end(), retval.begin() + ii * num_eq);
2021 }
2022 return retval;
2023 }
2024
2025
2026
2027
2028 template <typename TypeTag>
2029 int
2030 MultisegmentWell<TypeTag>::
2031 setPrimaryVars(std::vector<double>::const_iterator it)
2032 {
2033 const int num_seg = this->numberOfSegments();
2034 constexpr int num_eq = MSWEval::numWellEq;
2035 std::array<double, num_eq> tmp;
2036 for (int ii = 0; ii < num_seg; ++ii) {
2037 const auto start = it + num_seg * num_eq;
2038 std::copy(start, start + num_eq, tmp.begin());
2039 this->primary_variables_.setValue(ii, tmp);
2040 }
2041 return num_seg * num_eq;
2042 }
2043
2044} // namespace Opm
Definition AquiferInterface.hpp:35
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
Definition GroupState.hpp:34
Definition MultisegmentWell.hpp:37
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition MultisegmentWell_impl.hpp:220
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
Definition PerforationData.hpp:39