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