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