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