22#ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
23#define OPM_MULTISEGMENTWELL_IMPL_HEADER_INCLUDED
25#include <opm/simulators/wells/MultisegmentWell.hpp>
28#include <opm/common/Exceptions.hpp>
29#include <opm/common/OpmLog/OpmLog.hpp>
31#include <opm/input/eclipse/Schedule/MSW/Segment.hpp>
32#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
33#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
34#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
35#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
37#include <opm/input/eclipse/Units/Units.hpp>
39#include <opm/material/densead/EvaluationFormat.hpp>
41#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
42#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
43#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
49#if COMPILE_BDA_BRIDGE && (HAVE_CUDA || HAVE_OPENCL)
50#include <opm/simulators/linalg/bda/WellContributions.hpp>
57 template <
typename TypeTag>
58 MultisegmentWell<TypeTag>::
59 MultisegmentWell(
const Well& well,
60 const ParallelWellInfo<Scalar>& pw_info,
62 const ModelParameters& param,
63 const RateConverterType& rate_converter,
64 const int pvtRegionIdx,
65 const int num_components,
67 const int index_of_well,
68 const std::vector<PerforationData<Scalar>>& perf_data)
69 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
70 , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices>&>(*this))
72 , segment_fluid_initial_(this->numberOfSegments(), std::vector<Scalar>(this->num_components_, 0.0))
75 if constexpr (has_solvent) {
76 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
79 if constexpr (has_polymer) {
80 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
83 if constexpr (Base::has_energy) {
84 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
87 if constexpr (Base::has_foam) {
88 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
91 if constexpr (Base::has_brine) {
92 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
95 if constexpr (Base::has_watVapor) {
96 OPM_THROW(std::runtime_error,
"water evaporation is not supported by multisegment well yet");
99 if(this->rsRvInj() > 0) {
100 OPM_THROW(std::runtime_error,
101 "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
102 " \n See (WCONINJE item 10 / WCONHIST item 8)");
104 if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
105 OPM_THROW(std::runtime_error,
"water + gas case not supported by multisegment well yet");
108 this->thp_update_iterations =
true;
115 template <
typename TypeTag>
117 MultisegmentWell<TypeTag>::
118 init(
const PhaseUsage* phase_usage_arg,
119 const std::vector<Scalar>& depth_arg,
120 const Scalar gravity_arg,
121 const std::vector< Scalar >& B_avg,
122 const bool changed_to_open_this_step)
124 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
136 this->initMatrixAndVectors();
139 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
140 const int cell_idx = this->well_cells_[perf];
141 this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
149 template <
typename TypeTag>
151 MultisegmentWell<TypeTag>::
152 initPrimaryVariablesEvaluation()
154 this->primary_variables_.init();
161 template <
typename TypeTag>
163 MultisegmentWell<TypeTag>::
164 updatePrimaryVariables(
const Simulator& simulator,
165 const WellState<Scalar>& well_state,
166 DeferredLogger& deferred_logger)
168 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
169 this->primary_variables_.update(well_state, stop_or_zero_rate_target);
177 template <
typename TypeTag>
185 Base::updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
188 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
189 this->segments_.perforations(),
191 this->scaleSegmentPressuresWithBhp(well_state);
198 template <
typename TypeTag>
203 const std::vector<Scalar>& B_avg,
205 const bool relax_tolerance)
const
207 return this->MSWEval::getWellConvergence(well_state,
210 this->param_.max_residual_allowed_,
211 this->param_.tolerance_wells_,
212 this->param_.relaxed_tolerance_flow_well_,
213 this->param_.tolerance_pressure_ms_wells_,
214 this->param_.relaxed_tolerance_pressure_ms_well_,
216 this->wellIsStopped());
224 template <
typename TypeTag>
227 apply(
const BVector& x, BVector& Ax)
const
229 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
233 if (this->param_.matrix_add_well_contributions_) {
238 this->linSys_.
apply(x, Ax);
245 template <
typename TypeTag>
248 apply(BVector& r)
const
250 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
254 this->linSys_.
apply(r);
259 template <
typename TypeTag>
267 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
272 this->linSys_.recoverSolutionWell(x, xw);
273 updateWellState(simulator, xw, well_state, deferred_logger);
280 template <
typename TypeTag>
285 std::vector<Scalar>& well_potentials,
288 const auto [compute_potential, bhp_controlled_well] =
291 if (!compute_potential) {
295 debug_cost_counter_ = 0;
296 bool converged_implicit =
false;
297 if (this->param_.local_well_solver_control_switching_) {
298 converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
299 if (!converged_implicit) {
300 deferred_logger.debug(
"Implicit potential calculations failed for well "
301 + this->name() +
", reverting to original aproach.");
304 if (!converged_implicit) {
306 const auto& summaryState = simulator.vanguard().summaryState();
307 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
308 computeWellRatesAtBhpLimit(simulator, well_potentials, deferred_logger);
310 well_potentials = computeWellPotentialWithTHP(
311 well_state, simulator, deferred_logger);
314 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
315 + this->name() +
": " + std::to_string(debug_cost_counter_));
317 this->checkNegativeWellPotentials(well_potentials,
318 this->param_.check_well_operability_,
325 template<
typename TypeTag>
329 std::vector<Scalar>& well_flux,
332 if (this->well_ecl_.isInjector()) {
333 const auto controls = this->well_ecl_.injectionControls(simulator.
vanguard().summaryState());
334 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
336 const auto controls = this->well_ecl_.productionControls(simulator.
vanguard().summaryState());
337 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
341 template<
typename TypeTag>
343 MultisegmentWell<TypeTag>::
344 computeWellRatesWithBhp(
const Simulator& simulator,
346 std::vector<Scalar>& well_flux,
347 DeferredLogger& deferred_logger)
const
349 const int np = this->number_of_phases_;
351 well_flux.resize(np, 0.0);
352 const bool allow_cf = this->getAllowCrossFlow();
353 const int nseg = this->numberOfSegments();
354 const WellState<Scalar>& well_state = simulator.problem().wellModel().wellState();
355 const auto& ws = well_state.well(this->indexOfWell());
356 auto segments_copy = ws.segments;
357 segments_copy.scale_pressure(bhp);
358 const auto& segment_pressure = segments_copy.pressure;
359 for (
int seg = 0; seg < nseg; ++seg) {
360 for (
const int perf : this->segments_.perforations()[seg]) {
361 const int cell_idx = this->well_cells_[perf];
362 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
364 std::vector<Scalar> mob(this->num_components_, 0.);
365 getMobility(simulator, perf, mob, deferred_logger);
366 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
367 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
368 const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
369 const Scalar seg_pressure = segment_pressure[seg];
370 std::vector<Scalar> cq_s(this->num_components_, 0.);
371 Scalar perf_press = 0.0;
372 PerforationRates<Scalar> perf_rates;
373 computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
374 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
376 for(
int p = 0; p < np; ++p) {
377 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
381 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
385 template<
typename TypeTag>
387 MultisegmentWell<TypeTag>::
388 computeWellRatesWithBhpIterations(
const Simulator& simulator,
390 std::vector<Scalar>& well_flux,
391 DeferredLogger& deferred_logger)
const
395 MultisegmentWell<TypeTag> well_copy(*
this);
396 well_copy.resetDampening();
398 well_copy.debug_cost_counter_ = 0;
401 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
402 const auto& group_state = simulator.problem().wellModel().groupState();
403 auto& ws = well_state_copy.well(this->index_of_well_);
406 const auto& summary_state = simulator.vanguard().summaryState();
407 auto inj_controls = well_copy.well_ecl_.isInjector()
408 ? well_copy.well_ecl_.injectionControls(summary_state)
409 : Well::InjectionControls(0);
410 auto prod_controls = well_copy.well_ecl_.isProducer()
411 ? well_copy.well_ecl_.productionControls(summary_state) :
412 Well::ProductionControls(0);
415 if (well_copy.well_ecl_.isInjector()) {
416 inj_controls.bhp_limit = bhp;
417 ws.injection_cmode = Well::InjectorCMode::BHP;
419 prod_controls.bhp_limit = bhp;
420 ws.production_cmode = Well::ProducerCMode::BHP;
423 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
426 const int np = this->number_of_phases_;
428 for (
int phase = 0; phase < np; ++phase){
429 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
432 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
433 for (
int phase = 0; phase < np; ++phase) {
434 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
437 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
438 this->segments_.perforations(),
441 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
442 const double dt = simulator.timeStepSize();
444 well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
449 well_flux.resize(np, 0.0);
450 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
451 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
452 well_flux[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
454 debug_cost_counter_ += well_copy.debug_cost_counter_;
459 template<
typename TypeTag>
460 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
461 MultisegmentWell<TypeTag>::
462 computeWellPotentialWithTHP(
const WellState<Scalar>& well_state,
463 const Simulator& simulator,
464 DeferredLogger& deferred_logger)
const
466 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
467 const auto& summary_state = simulator.vanguard().summaryState();
469 const auto& well = this->well_ecl_;
470 if (well.isInjector()) {
471 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
472 if (bhp_at_thp_limit) {
473 const auto& controls = well.injectionControls(summary_state);
474 const Scalar bhp = std::min(*bhp_at_thp_limit,
475 static_cast<Scalar
>(controls.bhp_limit));
476 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
477 deferred_logger.debug(
"Converged thp based potential calculation for well "
478 + this->name() +
", at bhp = " + std::to_string(bhp));
480 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
481 "Failed in getting converged thp based potential calculation for well "
482 + this->name() +
". Instead the bhp based value is used");
483 const auto& controls = well.injectionControls(summary_state);
484 const Scalar bhp = controls.bhp_limit;
485 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
488 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
489 well_state, simulator, summary_state, deferred_logger);
490 if (bhp_at_thp_limit) {
491 const auto& controls = well.productionControls(summary_state);
492 const Scalar bhp = std::max(*bhp_at_thp_limit,
493 static_cast<Scalar
>(controls.bhp_limit));
494 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
495 deferred_logger.debug(
"Converged thp based potential calculation for well "
496 + this->name() +
", at bhp = " + std::to_string(bhp));
498 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
499 "Failed in getting converged thp based potential calculation for well "
500 + this->name() +
". Instead the bhp based value is used");
501 const auto& controls = well.productionControls(summary_state);
502 const Scalar bhp = controls.bhp_limit;
503 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
510 template<
typename TypeTag>
512 MultisegmentWell<TypeTag>::
513 computeWellPotentialsImplicit(
const Simulator& simulator,
514 std::vector<Scalar>& well_potentials,
515 DeferredLogger& deferred_logger)
const
520 MultisegmentWell<TypeTag> well_copy(*
this);
521 well_copy.debug_cost_counter_ = 0;
524 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
525 const auto& group_state = simulator.problem().wellModel().groupState();
526 auto& ws = well_state_copy.well(this->index_of_well_);
529 const auto& summary_state = simulator.vanguard().summaryState();
530 auto inj_controls = well_copy.well_ecl_.isInjector()
531 ? well_copy.well_ecl_.injectionControls(summary_state)
532 : Well::InjectionControls(0);
533 auto prod_controls = well_copy.well_ecl_.isProducer()
534 ? well_copy.well_ecl_.productionControls(summary_state)
535 : Well::ProductionControls(0);
538 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
540 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
543 const int np = this->number_of_phases_;
545 for (
int phase = 0; phase < np; ++phase){
546 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
549 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
550 for (
int phase = 0; phase < np; ++phase) {
551 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
554 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
555 this->segments_.perforations(),
558 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
559 const double dt = simulator.timeStepSize();
561 bool converged =
false;
562 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
563 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
565 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
569 well_potentials.clear();
570 well_potentials.resize(np, 0.0);
571 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
572 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
573 well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
575 debug_cost_counter_ += well_copy.debug_cost_counter_;
579 template <
typename TypeTag>
581 MultisegmentWell<TypeTag>::
582 solveEqAndUpdateWellState(
const Simulator& simulator,
583 WellState<Scalar>& well_state,
584 DeferredLogger& deferred_logger)
586 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
591 const BVectorWell dx_well = this->linSys_.solve();
593 updateWellState(simulator, dx_well, well_state, deferred_logger);
595 catch(
const NumericalProblem& exp) {
599 deferred_logger.problem(
"In MultisegmentWell::solveEqAndUpdateWellState for well "
600 + this->name() +
": "+exp.what());
609 template <
typename TypeTag>
611 MultisegmentWell<TypeTag>::
612 computePerfCellPressDiffs(
const Simulator& simulator)
614 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
616 std::vector<Scalar> kr(this->number_of_phases_, 0.0);
617 std::vector<Scalar> density(this->number_of_phases_, 0.0);
619 const int cell_idx = this->well_cells_[perf];
620 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
621 const auto& fs = intQuants.fluidState();
626 if (pu.phase_used[Water]) {
627 const int water_pos = pu.phase_pos[Water];
628 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
629 sum_kr += kr[water_pos];
630 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
633 if (pu.phase_used[Oil]) {
634 const int oil_pos = pu.phase_pos[Oil];
635 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
636 sum_kr += kr[oil_pos];
637 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
640 if (pu.phase_used[Gas]) {
641 const int gas_pos = pu.phase_pos[Gas];
642 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
643 sum_kr += kr[gas_pos];
644 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
647 assert(sum_kr != 0.);
650 Scalar average_density = 0.;
651 for (
int p = 0; p < this->number_of_phases_; ++p) {
652 average_density += kr[p] * density[p];
654 average_density /= sum_kr;
656 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
664 template <
typename TypeTag>
666 MultisegmentWell<TypeTag>::
667 computeInitialSegmentFluids(
const Simulator& simulator)
669 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
671 const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
672 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
673 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
682 template <
typename TypeTag>
684 MultisegmentWell<TypeTag>::
685 updateWellState(
const Simulator& simulator,
686 const BVectorWell& dwells,
687 WellState<Scalar>& well_state,
688 DeferredLogger& deferred_logger,
689 const Scalar relaxation_factor)
691 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
693 const Scalar dFLimit = this->param_.dwell_fraction_max_;
694 const Scalar max_pressure_change = this->param_.max_pressure_change_ms_wells_;
695 const bool stop_or_zero_rate_target =
696 this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
697 this->primary_variables_.updateNewton(dwells,
700 stop_or_zero_rate_target,
701 max_pressure_change);
703 const auto& summary_state = simulator.vanguard().summaryState();
704 this->primary_variables_.copyToWellState(*
this, getRefDensity(),
710 auto& ws = well_state.well(this->index_of_well_);
711 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
714 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
721 template <
typename TypeTag>
723 MultisegmentWell<TypeTag>::
724 calculateExplicitQuantities(
const Simulator& simulator,
725 const WellState<Scalar>& well_state,
726 DeferredLogger& deferred_logger)
728 updatePrimaryVariables(simulator, well_state, deferred_logger);
729 initPrimaryVariablesEvaluation();
730 computePerfCellPressDiffs(simulator);
731 computeInitialSegmentFluids(simulator);
738 template<
typename TypeTag>
740 MultisegmentWell<TypeTag>::
741 updateProductivityIndex(
const Simulator& simulator,
742 const WellProdIndexCalculator<Scalar>& wellPICalc,
743 WellState<Scalar>& well_state,
744 DeferredLogger& deferred_logger)
const
746 auto fluidState = [&simulator,
this](
const int perf)
748 const auto cell_idx = this->well_cells_[perf];
749 return simulator.model()
750 .intensiveQuantities(cell_idx, 0).fluidState();
753 const int np = this->number_of_phases_;
754 auto setToZero = [np](Scalar* x) ->
void
756 std::fill_n(x, np, 0.0);
759 auto addVector = [np](
const Scalar* src, Scalar* dest) ->
void
761 std::transform(src, src + np, dest, dest, std::plus<>{});
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();
771 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
772 auto subsetPerfID = 0;
774 for (
const auto& perf : *this->perf_data_){
775 auto allPerfID = perf.ecl_index;
777 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) -> Scalar
779 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
782 std::vector<Scalar> mob(this->num_components_, 0.0);
783 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
785 const auto& fs = fluidState(subsetPerfID);
788 if (this->isInjector()) {
789 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
790 mob, connPI, deferred_logger);
793 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
796 addVector(connPI, wellPI);
802 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
803 "Internal logic error in processing connections for PI/II");
810 template<
typename TypeTag>
811 typename MultisegmentWell<TypeTag>::Scalar
812 MultisegmentWell<TypeTag>::
813 connectionDensity(
const int globalConnIdx,
814 [[maybe_unused]]
const int openConnIdx)
const
819 const auto segNum = this->wellEcl()
820 .getConnections()[globalConnIdx].segment();
822 const auto segIdx = this->wellEcl()
823 .getSegments().segmentNumberToIndex(segNum);
825 return this->segments_.density(segIdx).value();
832 template<
typename TypeTag>
834 MultisegmentWell<TypeTag>::
835 addWellContributions(SparseMatrixAdapter& jacobian)
const
837 this->linSys_.extract(jacobian);
841 template<
typename TypeTag>
843 MultisegmentWell<TypeTag>::
844 addWellPressureEquations(PressureMatrix& jacobian,
845 const BVector& weights,
846 const int pressureVarIndex,
847 const bool use_well_weights,
848 const WellState<Scalar>& well_state)
const
851 this->linSys_.extractCPRPressureMatrix(jacobian,
861 template<
typename TypeTag>
862 template<
class Value>
864 MultisegmentWell<TypeTag>::
865 computePerfRate(
const Value& pressure_cell,
868 const std::vector<Value>& b_perfcells,
869 const std::vector<Value>& mob_perfcells,
870 const std::vector<Scalar>& Tw,
872 const Value& segment_pressure,
873 const Value& segment_density,
874 const bool& allow_cf,
875 const std::vector<Value>& cmix_s,
876 std::vector<Value>& cq_s,
878 PerforationRates<Scalar>& perf_rates,
879 DeferredLogger& deferred_logger)
const
882 const Value perf_seg_press_diff = this->gravity() * segment_density *
883 this->segments_.perforation_depth_diff(perf);
885 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
889 perf_press = segment_pressure + perf_seg_press_diff;
892 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
895 const Value drawdown = cell_press_at_perf - perf_press;
898 if (drawdown > 0.0) {
900 if (!allow_cf && this->isInjector()) {
905 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
906 const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
907 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
910 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
911 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
912 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
913 const Value cq_s_oil = cq_s[oilCompIdx];
914 const Value cq_s_gas = cq_s[gasCompIdx];
915 cq_s[gasCompIdx] += rs * cq_s_oil;
916 cq_s[oilCompIdx] += rv * cq_s_gas;
920 if (!allow_cf && this->isProducer()) {
925 Value total_mob = mob_perfcells[0];
926 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
927 total_mob += mob_perfcells[comp_idx];
931 Value volume_ratio = 0.0;
932 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
933 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
934 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
937 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
938 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
939 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
944 const Value d = 1.0 - rv * rs;
946 if (getValue(d) == 0.0) {
947 OPM_DEFLOG_PROBLEM(NumericalProblem,
948 fmt::format(
"Zero d value obtained for well {} "
949 "during flux calculation with rs {} and rv {}",
950 this->name(), rs, rv),
954 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
955 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
957 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
958 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
960 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
961 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
962 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
964 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
965 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
966 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
970 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
971 const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown);
972 Value cqt_is = cqt_i / volume_ratio;
973 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
978 if (this->isProducer()) {
979 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
980 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
981 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
990 const Scalar d = 1.0 - getValue(rv) * getValue(rs);
993 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
996 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
1001 template <
typename TypeTag>
1002 template<
class Value>
1004 MultisegmentWell<TypeTag>::
1005 computePerfRate(
const IntensiveQuantities& int_quants,
1006 const std::vector<Value>& mob_perfcells,
1007 const std::vector<Scalar>& Tw,
1010 const Value& segment_pressure,
1011 const bool& allow_cf,
1012 std::vector<Value>& cq_s,
1014 PerforationRates<Scalar>& perf_rates,
1015 DeferredLogger& deferred_logger)
const
1018 auto obtain = [
this](
const Eval& value)
1020 if constexpr (std::is_same_v<Value, Scalar>) {
1021 static_cast<void>(
this);
1022 return getValue(value);
1024 return this->extendEval(value);
1027 auto obtainN = [](
const auto& value)
1029 if constexpr (std::is_same_v<Value, Scalar>) {
1030 return getValue(value);
1035 const auto& fs = int_quants.fluidState();
1037 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
1038 const Value rs = obtain(fs.Rs());
1039 const Value rv = obtain(fs.Rv());
1042 std::vector<Value> b_perfcells(this->num_components_, 0.0);
1044 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1045 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1049 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1050 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
1053 std::vector<Value> cmix_s(this->numComponents(), 0.0);
1054 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
1055 cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
1058 this->computePerfRate(pressure_cell,
1066 obtainN(this->segments_.density(seg)),
1075 template <
typename TypeTag>
1077 MultisegmentWell<TypeTag>::
1078 computeSegmentFluidProperties(
const Simulator& simulator, DeferredLogger& deferred_logger)
1087 EvalWell temperature;
1088 EvalWell saltConcentration;
1094 int pvt_region_index;
1097 const int cell_idx = this->well_cells_[0];
1098 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1099 const auto& fs = intQuants.fluidState();
1100 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1101 saltConcentration = this->extendEval(fs.saltConcentration());
1102 pvt_region_index = fs.pvtRegionIndex();
1105 this->segments_.computeFluidProperties(temperature,
1107 this->primary_variables_,
1112 template <
typename TypeTag>
1113 template<
class Value>
1115 MultisegmentWell<TypeTag>::
1116 getMobility(
const Simulator& simulator,
1118 std::vector<Value>& mob,
1119 DeferredLogger& deferred_logger)
const
1121 auto obtain = [
this](
const Eval& value)
1123 if constexpr (std::is_same_v<Value, Scalar>) {
1124 static_cast<void>(
this);
1125 return getValue(value);
1127 return this->extendEval(value);
1131 WellInterface<TypeTag>::getMobility(simulator, perf, mob, obtain, deferred_logger);
1133 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
1134 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1135 const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
1136 const int seg = this->segmentNumberToIndex(con.segment());
1140 const Scalar segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
1141 const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
1142 * this->segments_.perforation_depth_diff(perf);
1143 const Scalar perf_press = segment_pres + perf_seg_press_diff;
1144 const Scalar multiplier = this->getInjMult(perf, segment_pres, perf_press, deferred_logger);
1145 for (std::size_t i = 0; i < mob.size(); ++i) {
1146 mob[i] *= multiplier;
1153 template<
typename TypeTag>
1154 typename MultisegmentWell<TypeTag>::Scalar
1155 MultisegmentWell<TypeTag>::
1156 getRefDensity()
const
1158 return this->segments_.getRefDensity();
1161 template<
typename TypeTag>
1163 MultisegmentWell<TypeTag>::
1164 checkOperabilityUnderBHPLimit(
const WellState<Scalar>& ,
1165 const Simulator& simulator,
1166 DeferredLogger& deferred_logger)
1168 const auto& summaryState = simulator.vanguard().summaryState();
1169 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1172 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1173 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1176 Scalar total_ipr_mass_rate = 0.0;
1177 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1179 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1183 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1184 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1186 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1187 total_ipr_mass_rate += ipr_rate * rho;
1189 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1190 this->operability_status_.operable_under_only_bhp_limit =
false;
1194 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1198 std::vector<Scalar> well_rates_bhp_limit;
1199 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1201 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1202 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1204 this->getRefDensity(),
1205 this->wellEcl().alq_value(summaryState),
1208 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1209 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1220 this->operability_status_.operable_under_only_bhp_limit =
true;
1221 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1227 template<
typename TypeTag>
1229 MultisegmentWell<TypeTag>::
1230 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
1235 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1236 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1238 const int nseg = this->numberOfSegments();
1239 std::vector<Scalar> seg_dp(nseg, 0.0);
1240 for (
int seg = 0; seg < nseg; ++seg) {
1242 const Scalar dp = this->getSegmentDp(seg,
1243 this->segments_.density(seg).value(),
1246 for (
const int perf : this->segments_.perforations()[seg]) {
1247 std::vector<Scalar> mob(this->num_components_, 0.0);
1250 getMobility(simulator, perf, mob, deferred_logger);
1252 const int cell_idx = this->well_cells_[perf];
1253 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
1254 const auto& fs = int_quantities.fluidState();
1256 const Scalar perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1258 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1259 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1262 std::vector<Scalar> b_perf(this->num_components_);
1263 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1264 if (!FluidSystem::phaseIsActive(phase)) {
1267 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1268 b_perf[comp_idx] = fs.invB(phase).value();
1272 const Scalar h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1273 const Scalar pressure_diff = pressure_cell - h_perf;
1276 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1277 deferred_logger.debug(
"CROSSFLOW_IPR",
1278 "cross flow found when updateIPR for well " + this->name());
1282 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
1283 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1284 const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
1285 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
1286 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
1287 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1288 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
1289 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1290 ipr_b_perf[comp_idx] += tw_mob;
1294 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1295 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1296 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1297 const Scalar rs = (fs.Rs()).value();
1298 const Scalar rv = (fs.Rv()).value();
1300 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1301 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1303 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1304 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1306 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1307 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1309 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1310 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1313 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1314 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1315 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1321 template<
typename TypeTag>
1323 MultisegmentWell<TypeTag>::
1324 updateIPRImplicit(
const Simulator& simulator,
1325 WellState<Scalar>& well_state,
1326 DeferredLogger& deferred_logger)
1335 auto rates = well_state.well(this->index_of_well_).surface_rates;
1337 for (std::size_t p = 0; p < rates.size(); ++p) {
1338 zero_rates &= rates[p] == 0.0;
1340 auto& ws = well_state.well(this->index_of_well_);
1342 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
1343 deferred_logger.debug(msg);
1355 const auto& group_state = simulator.problem().wellModel().groupState();
1357 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
1358 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
1360 auto inj_controls = Well::InjectionControls(0);
1361 auto prod_controls = Well::ProductionControls(0);
1362 prod_controls.addControl(Well::ProducerCMode::BHP);
1363 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
1366 const auto cmode = ws.production_cmode;
1367 ws.production_cmode = Well::ProducerCMode::BHP;
1368 const double dt = simulator.timeStepSize();
1369 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1371 BVectorWell rhs(this->numberOfSegments());
1373 rhs[0][SPres] = -1.0;
1375 const BVectorWell x_well = this->linSys_.solve(rhs);
1376 constexpr int num_eq = MSWEval::numWellEq;
1377 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
1378 const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
1379 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
1380 for (
size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
1382 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
1384 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
1387 ws.production_cmode = cmode;
1390 template<
typename TypeTag>
1392 MultisegmentWell<TypeTag>::
1393 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
1394 const WellState<Scalar>& well_state,
1395 DeferredLogger& deferred_logger)
1397 const auto& summaryState = simulator.vanguard().summaryState();
1398 const auto obtain_bhp = this->isProducer()
1399 ? computeBhpAtThpLimitProd(
1400 well_state, simulator, summaryState, deferred_logger)
1401 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1404 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1406 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1407 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1409 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1410 if (this->isProducer() && *obtain_bhp < thp_limit) {
1411 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1412 +
" bars is SMALLER than thp limit "
1413 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1414 +
" bars as a producer for well " + this->name();
1415 deferred_logger.debug(msg);
1417 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1418 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1419 +
" bars is LARGER than thp limit "
1420 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1421 +
" bars as a injector for well " + this->name();
1422 deferred_logger.debug(msg);
1427 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1428 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1429 if (!this->wellIsStopped()) {
1430 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1431 deferred_logger.debug(
" could not find bhp value at thp limit "
1432 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1433 +
" bar for well " + this->name() +
", the well might need to be closed ");
1442 template<
typename TypeTag>
1444 MultisegmentWell<TypeTag>::
1445 iterateWellEqWithControl(
const Simulator& simulator,
1447 const Well::InjectionControls& inj_controls,
1448 const Well::ProductionControls& prod_controls,
1449 WellState<Scalar>& well_state,
1450 const GroupState<Scalar>& group_state,
1451 DeferredLogger& deferred_logger)
1453 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1455 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1459 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1464 std::vector<std::vector<Scalar> > residual_history;
1465 std::vector<Scalar> measure_history;
1468 Scalar relaxation_factor = 1.;
1469 const Scalar min_relaxation_factor = 0.6;
1470 bool converged =
false;
1471 int stagnate_count = 0;
1472 bool relax_convergence =
false;
1473 this->regularize_ =
false;
1474 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1476 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1478 BVectorWell dx_well;
1480 dx_well = this->linSys_.solve();
1482 catch(
const NumericalProblem& exp) {
1486 deferred_logger.problem(
"In MultisegmentWell::iterateWellEqWithControl for well "
1487 + this->name() +
": "+exp.what());
1491 if (it > this->param_.strict_inner_iter_wells_) {
1492 relax_convergence =
true;
1493 this->regularize_ =
true;
1496 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1497 if (report.converged()) {
1504 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1508 residual_history.push_back(residuals);
1509 measure_history.push_back(this->getResidualMeasureValue(well_state,
1510 residual_history[it],
1511 this->param_.tolerance_wells_,
1512 this->param_.tolerance_pressure_ms_wells_,
1517 bool is_oscillate =
false;
1518 bool is_stagnate =
false;
1520 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1524 if (is_oscillate || is_stagnate) {
1526 std::ostringstream sstr;
1527 if (relaxation_factor == min_relaxation_factor) {
1530 if (stagnate_count == 6) {
1531 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1532 const auto reportStag = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger,
true);
1533 if (reportStag.converged()) {
1535 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1536 deferred_logger.debug(sstr.str());
1543 const Scalar reduction_mutliplier = 0.9;
1544 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1548 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1552 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1554 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1556 this->regularize_ =
true;
1557 deferred_logger.debug(sstr.str());
1559 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1560 initPrimaryVariablesEvaluation();
1565 std::ostringstream sstr;
1566 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1567 if (relax_convergence)
1568 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1569 deferred_logger.debug(sstr.str());
1571 std::ostringstream sstr;
1572 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1573#define EXTRA_DEBUG_MSW 0
1575 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1576 for (
int i = 0; i < it; ++i) {
1577 const auto& residual = residual_history[i];
1578 sstr <<
" residual at " << i <<
"th iteration ";
1579 for (
const auto& res : residual) {
1582 sstr <<
" " << measure_history[i] <<
" \n";
1585#undef EXTRA_DEBUG_MSW
1586 deferred_logger.debug(sstr.str());
1593 template<
typename TypeTag>
1595 MultisegmentWell<TypeTag>::
1596 iterateWellEqWithSwitching(
const Simulator& simulator,
1598 const Well::InjectionControls& inj_controls,
1599 const Well::ProductionControls& prod_controls,
1600 WellState<Scalar>& well_state,
1601 const GroupState<Scalar>& group_state,
1602 DeferredLogger& deferred_logger,
1603 const bool fixed_control ,
1604 const bool fixed_status )
1606 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1610 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1615 std::vector<std::vector<Scalar> > residual_history;
1616 std::vector<Scalar> measure_history;
1619 Scalar relaxation_factor = 1.;
1620 const Scalar min_relaxation_factor = 0.6;
1621 bool converged =
false;
1622 [[maybe_unused]]
int stagnate_count = 0;
1623 bool relax_convergence =
false;
1624 this->regularize_ =
false;
1625 const auto& summary_state = simulator.vanguard().summaryState();
1630 const int min_its_after_switch = 3;
1631 int its_since_last_switch = min_its_after_switch;
1632 int switch_count= 0;
1634 const auto well_status_orig = this->wellStatus_;
1635 const auto operability_orig = this->operability_status_;
1636 auto well_status_cur = well_status_orig;
1638 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
1639 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
1641 const bool allow_switching = !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
1642 (!fixed_control || !fixed_status) && allow_open;
1643 bool changed =
false;
1644 bool final_check =
false;
1646 this->operability_status_.resetOperability();
1647 this->operability_status_.solvable =
true;
1649 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1650 its_since_last_switch++;
1651 if (allow_switching && its_since_last_switch >= min_its_after_switch){
1652 const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
1653 changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
1654 inj_controls, prod_controls, wqTotal,
1655 deferred_logger, fixed_control, fixed_status);
1657 its_since_last_switch = 0;
1659 if (well_status_cur != this->wellStatus_) {
1660 well_status_cur = this->wellStatus_;
1663 if (!changed && final_check) {
1666 final_check =
false;
1670 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1672 const BVectorWell dx_well = this->linSys_.solve();
1674 if (it > this->param_.strict_inner_iter_wells_) {
1675 relax_convergence =
true;
1676 this->regularize_ =
true;
1679 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1680 converged = report.converged();
1684 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1686 its_since_last_switch = min_its_after_switch;
1694 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1698 residual_history.push_back(residuals);
1702 measure_history.push_back(this->getResidualMeasureValue(well_state,
1703 residual_history[it],
1704 this->param_.tolerance_wells_,
1705 this->param_.tolerance_pressure_ms_wells_,
1708 bool is_oscillate =
false;
1709 bool is_stagnate =
false;
1711 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1715 if (is_oscillate || is_stagnate) {
1717 std::string message;
1718 if (relaxation_factor == min_relaxation_factor) {
1721 fmt::format_to(std::back_inserter(message),
" Well {} observes severe stagnation and/or oscillation."
1722 " We relax the tolerance and check for convergence. \n", this->name());
1723 const auto reportStag = getWellConvergence(simulator, well_state, Base::B_avg_,
1724 deferred_logger,
true);
1725 if (reportStag.converged()) {
1727 fmt::format_to(std::back_inserter(message),
" Well {} manages to get converged with relaxed tolerances in {} inner iterations", this->name(), it);
1728 deferred_logger.debug(message);
1735 constexpr Scalar reduction_mutliplier = 0.9;
1736 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1740 fmt::format_to(std::back_inserter(message),
" well {} observes stagnation in inner iteration {}\n", this->name(), it);
1743 fmt::format_to(std::back_inserter(message),
" well {} observes oscillation in inner iteration {}\n", this->name(), it);
1745 fmt::format_to(std::back_inserter(message),
" relaxation_factor is {} now\n", relaxation_factor);
1747 this->regularize_ =
true;
1748 deferred_logger.debug(message);
1751 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1752 initPrimaryVariablesEvaluation();
1756 if (allow_switching){
1758 const bool is_stopped = this->wellIsStopped();
1759 if (this->wellHasTHPConstraints(summary_state)){
1760 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
1761 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
1763 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
1766 std::string message = fmt::format(
" Well {} converged in {} inner iterations ("
1767 "{} control/status switches).", this->name(), it, switch_count);
1768 if (relax_convergence) {
1769 message.append(fmt::format(
" (A relaxed tolerance was used after {} iterations)",
1770 this->param_.strict_inner_iter_wells_));
1772 deferred_logger.debug(message);
1774 this->wellStatus_ = well_status_orig;
1775 this->operability_status_ = operability_orig;
1776 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
1777 "{} control/status switches).", this->name(), it, switch_count);
1778 deferred_logger.debug(message);
1779 this->primary_variables_.outputLowLimitPressureSegments(deferred_logger);
1786 template<
typename TypeTag>
1788 MultisegmentWell<TypeTag>::
1789 assembleWellEqWithoutIteration(
const Simulator& simulator,
1791 const Well::InjectionControls& inj_controls,
1792 const Well::ProductionControls& prod_controls,
1793 WellState<Scalar>& well_state,
1794 const GroupState<Scalar>& group_state,
1795 DeferredLogger& deferred_logger)
1797 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1800 this->segments_.updateUpwindingSegments(this->primary_variables_);
1803 computeSegmentFluidProperties(simulator, deferred_logger);
1806 this->linSys_.clear();
1808 auto& ws = well_state.well(this->index_of_well_);
1809 ws.phase_mixing_rates.fill(0.0);
1816 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1818 const int nseg = this->numberOfSegments();
1820 for (
int seg = 0; seg < nseg; ++seg) {
1824 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg);
1829 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1831 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1832 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1833 - segment_fluid_initial_[seg][comp_idx]) / dt;
1834 MultisegmentWellAssemble(*this).
1835 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1840 const int seg_upwind = this->segments_.upwinding_segment(seg);
1841 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1842 const EvalWell segment_rate =
1843 this->primary_variables_.getSegmentRateUpwinding(seg,
1846 this->well_efficiency_factor_;
1847 MultisegmentWellAssemble(*this).
1848 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1854 for (
const int inlet : this->segments_.inlets()[seg]) {
1855 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1856 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1857 const EvalWell inlet_rate =
1858 this->primary_variables_.getSegmentRateUpwinding(inlet,
1861 this->well_efficiency_factor_;
1862 MultisegmentWellAssemble(*this).
1863 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1869 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1870 auto& perf_data = ws.perf_data;
1871 auto& perf_rates = perf_data.phase_rates;
1872 auto& perf_press_state = perf_data.pressure;
1873 for (
const int perf : this->segments_.perforations()[seg]) {
1874 const int cell_idx = this->well_cells_[perf];
1875 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
1876 std::vector<EvalWell> mob(this->num_components_, 0.0);
1877 getMobility(simulator, perf, mob, deferred_logger);
1878 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
1879 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1880 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
1881 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1882 EvalWell perf_press;
1883 PerforationRates<Scalar> perfRates;
1884 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
1885 allow_cf, cq_s, perf_press, perfRates, deferred_logger);
1888 if (this->isProducer()) {
1889 ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas;
1890 ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil;
1891 perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perfRates.dis_gas;
1892 perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perfRates.vap_oil;
1896 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1897 perf_rates[perf*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1899 perf_press_state[perf] = perf_press.value();
1901 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1903 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1905 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1907 MultisegmentWellAssemble(*this).
1908 assemblePerforationEq(seg, perf, comp_idx, cq_s_effective, this->linSys_);
1914 const auto& summaryState = simulator.vanguard().summaryState();
1915 const Schedule& schedule = simulator.vanguard().schedule();
1916 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1917 MultisegmentWellAssemble(*this).
1918 assembleControlEq(well_state,
1925 this->primary_variables_,
1927 stopped_or_zero_target,
1930 const UnitSystem& unit_system = simulator.vanguard().eclState().getDeckUnitSystem();
1931 const auto& summary_state = simulator.vanguard().summaryState();
1932 this->assemblePressureEq(seg, unit_system, well_state, summary_state, this->param_.use_average_density_ms_wells_, deferred_logger);
1936 this->linSys_.createSolver();
1942 template<
typename TypeTag>
1944 MultisegmentWell<TypeTag>::
1945 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1947 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1951 template<
typename TypeTag>
1953 MultisegmentWell<TypeTag>::
1954 allDrawDownWrongDirection(
const Simulator& simulator)
const
1956 bool all_drawdown_wrong_direction =
true;
1957 const int nseg = this->numberOfSegments();
1959 for (
int seg = 0; seg < nseg; ++seg) {
1960 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1961 for (
const int perf : this->segments_.perforations()[seg]) {
1963 const int cell_idx = this->well_cells_[perf];
1964 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1965 const auto& fs = intQuants.fluidState();
1968 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1970 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1972 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1973 const Scalar perf_press = pressure_cell - cell_perf_press_diff;
1976 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1981 if ( (drawdown < 0. && this->isInjector()) ||
1982 (drawdown > 0. && this->isProducer()) ) {
1983 all_drawdown_wrong_direction =
false;
1989 return all_drawdown_wrong_direction;
1995 template<
typename TypeTag>
1997 MultisegmentWell<TypeTag>::
1998 updateWaterThroughput(
const double , WellState<Scalar>& )
const
2006 template<
typename TypeTag>
2007 typename MultisegmentWell<TypeTag>::EvalWell
2008 MultisegmentWell<TypeTag>::
2009 getSegmentSurfaceVolume(
const Simulator& simulator,
const int seg_idx)
const
2011 EvalWell temperature;
2012 EvalWell saltConcentration;
2013 int pvt_region_index;
2017 const int cell_idx = this->well_cells_[0];
2018 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2019 const auto& fs = intQuants.fluidState();
2020 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
2021 saltConcentration = this->extendEval(fs.saltConcentration());
2022 pvt_region_index = fs.pvtRegionIndex();
2025 return this->segments_.getSurfaceVolume(temperature,
2027 this->primary_variables_,
2033 template<
typename TypeTag>
2034 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2035 MultisegmentWell<TypeTag>::
2036 computeBhpAtThpLimitProd(
const WellState<Scalar>& well_state,
2037 const Simulator& simulator,
2038 const SummaryState& summary_state,
2039 DeferredLogger& deferred_logger)
const
2041 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
2044 this->getALQ(well_state),
2050 template<
typename TypeTag>
2051 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2052 MultisegmentWell<TypeTag>::
2053 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2054 const SummaryState& summary_state,
2055 const Scalar alq_value,
2056 DeferredLogger& deferred_logger)
const
2059 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2065 std::vector<Scalar> rates(3);
2066 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2070 auto bhpAtLimit = WellBhpThpCalculator(*this).
2071 computeBhpAtThpLimitProd(frates,
2073 this->maxPerfPress(simulator),
2074 this->getRefDensity(),
2076 this->getTHPConstraint(summary_state),
2082 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2086 std::vector<Scalar> rates(3);
2087 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2091 return WellBhpThpCalculator(*this).
2092 computeBhpAtThpLimitProd(fratesIter,
2094 this->maxPerfPress(simulator),
2095 this->getRefDensity(),
2097 this->getTHPConstraint(summary_state),
2101 template<
typename TypeTag>
2102 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2103 MultisegmentWell<TypeTag>::
2104 computeBhpAtThpLimitInj(
const Simulator& simulator,
2105 const SummaryState& summary_state,
2106 DeferredLogger& deferred_logger)
const
2109 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2115 std::vector<Scalar> rates(3);
2116 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2120 auto bhpAtLimit = WellBhpThpCalculator(*this).
2121 computeBhpAtThpLimitInj(frates,
2123 this->getRefDensity(),
2132 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2136 std::vector<Scalar> rates(3);
2137 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2141 return WellBhpThpCalculator(*this).
2142 computeBhpAtThpLimitInj(fratesIter,
2144 this->getRefDensity(),
2155 template<
typename TypeTag>
2156 typename MultisegmentWell<TypeTag>::Scalar
2157 MultisegmentWell<TypeTag>::
2158 maxPerfPress(
const Simulator& simulator)
const
2160 Scalar max_pressure = 0.0;
2161 const int nseg = this->numberOfSegments();
2162 for (
int seg = 0; seg < nseg; ++seg) {
2163 for (
const int perf : this->segments_.perforations()[seg]) {
2164 const int cell_idx = this->well_cells_[perf];
2165 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2166 const auto& fs = int_quants.fluidState();
2167 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2168 max_pressure = std::max(max_pressure, pressure_cell);
2171 return max_pressure;
2178 template<
typename TypeTag>
2179 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2185 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
2186 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2187 const int nseg = this->numberOfSegments();
2188 for (
int seg = 0; seg < nseg; ++seg) {
2190 const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
2191 for (
const int perf : this->segments_.perforations()[seg]) {
2192 const int cell_idx = this->well_cells_[perf];
2193 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2194 std::vector<Scalar> mob(this->num_components_, 0.0);
2195 getMobility(simulator, perf, mob, deferred_logger);
2196 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
2197 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2198 const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
2199 std::vector<Scalar> cq_s(this->num_components_, 0.0);
2200 Scalar perf_press = 0.0;
2202 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
2203 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
2204 for (
int comp = 0; comp < this->num_components_; ++comp) {
2205 well_q_s[comp] += cq_s[comp];
2213 template <
typename TypeTag>
2214 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2218 const int num_seg = this->numberOfSegments();
2219 constexpr int num_eq = MSWEval::numWellEq;
2220 std::vector<Scalar> retval(num_seg * num_eq);
2221 for (
int ii = 0; ii < num_seg; ++ii) {
2222 const auto& pv = this->primary_variables_.value(ii);
2223 std::copy(pv.begin(), pv.end(), retval.begin() + ii * num_eq);
2231 template <
typename TypeTag>
2233 MultisegmentWell<TypeTag>::
2234 setPrimaryVars(
typename std::vector<Scalar>::const_iterator it)
2236 const int num_seg = this->numberOfSegments();
2237 constexpr int num_eq = MSWEval::numWellEq;
2238 std::array<Scalar, num_eq> tmp;
2239 for (
int ii = 0; ii < num_seg; ++ii) {
2240 const auto start = it + ii * num_eq;
2241 std::copy(start, start + num_eq, tmp.begin());
2242 this->primary_variables_.setValue(ii, tmp);
2244 return num_seg * num_eq;
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 WellInterfaceFluidSystem.hpp:43
Definition MultisegmentWell.hpp:36
void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition MultisegmentWell_impl.hpp:227
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:260
Definition WellTest.hpp:37
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:62
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
Definition PerforationData.hpp:41