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