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