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