21#include <opm/common/Exceptions.hpp>
22#include <opm/common/OpmLog/OpmLog.hpp>
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>
30#include <opm/input/eclipse/Units/Units.hpp>
32#include <opm/material/densead/EvaluationFormat.hpp>
34#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
35#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
36#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
42#if HAVE_CUDA || HAVE_OPENCL
43#include <opm/simulators/linalg/bda/WellContributions.hpp>
50 template <
typename TypeTag>
51 MultisegmentWell<TypeTag>::
52 MultisegmentWell(
const Well& well,
53 const ParallelWellInfo& pw_info,
55 const ModelParameters& param,
56 const RateConverterType& rate_converter,
57 const int pvtRegionIdx,
58 const int num_components,
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))
65 , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
68 if constexpr (has_solvent) {
69 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
72 if constexpr (has_polymer) {
73 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
76 if constexpr (Base::has_energy) {
77 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
80 if constexpr (Base::has_foam) {
81 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
84 if constexpr (Base::has_brine) {
85 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
88 if constexpr (Base::has_watVapor) {
89 OPM_THROW(std::runtime_error,
"water evaporation is not supported by multisegment well yet");
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)");
97 if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
98 OPM_THROW(std::runtime_error,
"water + gas case not supported by multisegment well yet");
101 this->thp_update_iterations =
true;
108 template <
typename TypeTag>
110 MultisegmentWell<TypeTag>::
111 init(
const PhaseUsage* phase_usage_arg,
112 const std::vector<double>& depth_arg,
113 const double gravity_arg,
115 const std::vector< Scalar >& B_avg,
116 const bool changed_to_open_this_step)
118 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
130 this->initMatrixAndVectors(num_cells);
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];
143 template <
typename TypeTag>
145 MultisegmentWell<TypeTag>::
146 initPrimaryVariablesEvaluation()
148 this->primary_variables_.init();
155 template <
typename TypeTag>
157 MultisegmentWell<TypeTag>::
158 updatePrimaryVariables(
const SummaryState& summary_state,
159 const WellState& well_state,
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);
171 template <
typename TypeTag>
179 Base::updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
182 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
183 this->segments_.perforations(),
185 this->scaleSegmentPressuresWithBhp(well_state);
192 template <
typename TypeTag>
197 const std::vector<double>& B_avg,
199 const bool relax_tolerance)
const
201 return this->MSWEval::getWellConvergence(well_state,
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_,
210 this->wellIsStopped());
218 template <
typename TypeTag>
221 apply(
const BVector& x, BVector& Ax)
const
223 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
227 if (this->param_.matrix_add_well_contributions_) {
232 this->linSys_.
apply(x, Ax);
239 template <
typename TypeTag>
242 apply(BVector& r)
const
244 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
248 this->linSys_.
apply(r);
253 template <
typename TypeTag>
261 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
266 this->linSys_.recoverSolutionWell(x, xw);
267 updateWellState(summary_state, xw, well_state, deferred_logger);
274 template <
typename TypeTag>
279 std::vector<double>& well_potentials,
282 const auto [compute_potential, bhp_controlled_well] =
283 this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
285 if (!compute_potential) {
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.");
298 if (!converged_implicit) {
300 const auto& summaryState = simulator.vanguard().summaryState();
301 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
302 computeWellRatesAtBhpLimit(simulator, well_potentials, deferred_logger);
304 well_potentials = computeWellPotentialWithTHP(
305 well_state, simulator, deferred_logger);
308 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
309 + this->name() +
": " + std::to_string(debug_cost_counter_));
311 this->checkNegativeWellPotentials(well_potentials,
312 this->param_.check_well_operability_,
319 template<
typename TypeTag>
323 std::vector<double>& well_flux,
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);
330 const auto controls = this->well_ecl_.productionControls(simulator.vanguard().summaryState());
331 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
335 template<
typename TypeTag>
337 MultisegmentWell<TypeTag>::
338 computeWellRatesWithBhp(
const Simulator& simulator,
340 std::vector<double>& well_flux,
341 DeferredLogger& deferred_logger)
const
344 const int np = this->number_of_phases_;
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, 0);
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);
371 for(
int p = 0; p < np; ++p) {
372 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
376 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
380 template<
typename TypeTag>
382 MultisegmentWell<TypeTag>::
383 computeWellRatesWithBhpIterations(
const Simulator& simulator,
385 std::vector<double>& well_flux,
386 DeferredLogger& deferred_logger)
const
390 MultisegmentWell<TypeTag> well_copy(*
this);
391 well_copy.debug_cost_counter_ = 0;
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_);
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);
408 if (well_copy.well_ecl_.isInjector()) {
409 inj_controls.bhp_limit = bhp;
410 ws.injection_cmode = Well::InjectorCMode::BHP;
412 prod_controls.bhp_limit = bhp;
413 ws.production_cmode = Well::ProducerCMode::BHP;
416 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
419 const int np = this->number_of_phases_;
421 for (
int phase = 0; phase < np; ++phase){
422 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
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];
430 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
431 this->segments_.perforations(),
434 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
435 const double dt = simulator.timeStepSize();
437 well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
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();
447 debug_cost_counter_ += well_copy.debug_cost_counter_;
452 template<
typename TypeTag>
454 MultisegmentWell<TypeTag>::
455 computeWellPotentialWithTHP(
456 const WellState& well_state,
457 const Simulator& simulator,
458 DeferredLogger& deferred_logger)
const
460 std::vector<double> potentials(this->number_of_phases_, 0.0);
461 const auto& summary_state = simulator.vanguard().summaryState();
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));
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);
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));
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);
502 template<
typename TypeTag>
504 MultisegmentWell<TypeTag>::
505 computeWellPotentialsImplicit(
const Simulator& simulator,
506 std::vector<double>& well_potentials,
507 DeferredLogger& deferred_logger)
const
512 MultisegmentWell<TypeTag> well_copy(*
this);
513 well_copy.debug_cost_counter_ = 0;
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_);
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);
530 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
532 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
535 const int np = this->number_of_phases_;
537 for (
int phase = 0; phase < np; ++phase){
538 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
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];
546 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
547 this->segments_.perforations(),
550 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
551 const double dt = simulator.timeStepSize();
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);
557 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
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();
567 debug_cost_counter_ += well_copy.debug_cost_counter_;
571 template <
typename TypeTag>
573 MultisegmentWell<TypeTag>::
574 solveEqAndUpdateWellState(
const SummaryState& summary_state,
575 WellState& well_state,
576 DeferredLogger& deferred_logger)
578 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
583 const BVectorWell dx_well = this->linSys_.solve();
585 updateWellState(summary_state, dx_well, well_state, deferred_logger);
587 catch(
const NumericalProblem& exp) {
591 deferred_logger.problem(
"In MultisegmentWell::solveEqAndUpdateWellState for well "
592 + this->name() +
": "+exp.what());
601 template <
typename TypeTag>
603 MultisegmentWell<TypeTag>::
604 computePerfCellPressDiffs(
const Simulator& simulator)
606 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
608 std::vector<double> kr(this->number_of_phases_, 0.0);
609 std::vector<double> density(this->number_of_phases_, 0.0);
611 const int cell_idx = this->well_cells_[perf];
612 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
613 const auto& fs = intQuants.fluidState();
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();
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();
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();
639 assert(sum_kr != 0.);
642 double average_density = 0.;
643 for (
int p = 0; p < this->number_of_phases_; ++p) {
644 average_density += kr[p] * density[p];
646 average_density /= sum_kr;
648 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
656 template <
typename TypeTag>
658 MultisegmentWell<TypeTag>::
659 computeInitialSegmentFluids(
const Simulator& simulator)
661 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
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();
674 template <
typename TypeTag>
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)
683 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
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,
691 stop_or_zero_rate_target,
692 max_pressure_change);
694 this->primary_variables_.copyToWellState(*
this, getRefDensity(), stop_or_zero_rate_target,
695 well_state, summary_state, deferred_logger);
698 auto& ws = well_state.well(this->index_of_well_);
699 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
702 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
709 template <
typename TypeTag>
711 MultisegmentWell<TypeTag>::
712 calculateExplicitQuantities(
const Simulator& simulator,
713 const WellState& well_state,
714 DeferredLogger& deferred_logger)
716 const auto& summary_state = simulator.vanguard().summaryState();
717 updatePrimaryVariables(summary_state, well_state, deferred_logger);
718 initPrimaryVariablesEvaluation();
719 computePerfCellPressDiffs(simulator);
720 computeInitialSegmentFluids(simulator);
727 template<
typename TypeTag>
729 MultisegmentWell<TypeTag>::
730 updateProductivityIndex(
const Simulator& simulator,
731 const WellProdIndexCalculator& wellPICalc,
732 WellState& well_state,
733 DeferredLogger& deferred_logger)
const
735 auto fluidState = [&simulator,
this](
const int perf)
737 const auto cell_idx = this->well_cells_[perf];
738 return simulator.model()
739 .intensiveQuantities(cell_idx, 0).fluidState();
742 const int np = this->number_of_phases_;
743 auto setToZero = [np](
double* x) ->
void
745 std::fill_n(x, np, 0.0);
748 auto addVector = [np](
const double* src,
double* dest) ->
void
750 std::transform(src, src + np, dest, dest, std::plus<>{});
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();
760 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
761 auto subsetPerfID = 0;
763 for (
const auto& perf : *this->perf_data_){
764 auto allPerfID = perf.ecl_index;
766 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
768 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
771 std::vector<Scalar> mob(this->num_components_, 0.0);
772 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
774 const auto& fs = fluidState(subsetPerfID);
777 if (this->isInjector()) {
778 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
779 mob, connPI, deferred_logger);
782 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
785 addVector(connPI, wellPI);
791 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
792 "Internal logic error in processing connections for PI/II");
799 template<
typename TypeTag>
801 MultisegmentWell<TypeTag>::
802 connectionDensity(
const int globalConnIdx,
803 [[maybe_unused]]
const int openConnIdx)
const
808 const auto segNum = this->wellEcl()
809 .getConnections()[globalConnIdx].segment();
811 const auto segIdx = this->wellEcl()
812 .getSegments().segmentNumberToIndex(segNum);
814 return this->segments_.density(segIdx).value();
821 template<
typename TypeTag>
823 MultisegmentWell<TypeTag>::
824 addWellContributions(SparseMatrixAdapter& jacobian)
const
826 this->linSys_.extract(jacobian);
830 template<
typename TypeTag>
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
840 this->linSys_.extractCPRPressureMatrix(jacobian,
850 template<
typename TypeTag>
851 template<
class Value>
853 MultisegmentWell<TypeTag>::
854 computePerfRate(
const Value& pressure_cell,
857 const std::vector<Value>& b_perfcells,
858 const std::vector<Value>& mob_perfcells,
859 const std::vector<Scalar>& Tw,
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,
867 PerforationRates& perf_rates,
868 DeferredLogger& deferred_logger)
const
871 const Value perf_seg_press_diff = this->gravity() * segment_density *
872 this->segments_.perforation_depth_diff(perf);
874 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
878 perf_press = segment_pressure + perf_seg_press_diff;
881 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
884 const Value drawdown = cell_press_at_perf - perf_press;
887 if (drawdown > 0.0) {
889 if (!allow_cf && this->isInjector()) {
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;
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;
909 if (!allow_cf && this->isProducer()) {
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];
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];
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);
933 const Value d = 1.0 - rv * rs;
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),
943 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
944 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
946 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
947 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
949 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
950 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
951 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
953 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
954 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
955 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
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;
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);
979 const double d = 1.0 - getValue(rv) * getValue(rs);
982 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
985 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
990 template <
typename TypeTag>
991 template<
class Value>
993 MultisegmentWell<TypeTag>::
994 computePerfRate(
const IntensiveQuantities& int_quants,
995 const std::vector<Value>& mob_perfcells,
996 const std::vector<Scalar>& Tw,
999 const Value& segment_pressure,
1000 const bool& allow_cf,
1001 std::vector<Value>& cq_s,
1003 PerforationRates& perf_rates,
1004 DeferredLogger& deferred_logger)
const
1007 auto obtain = [
this](
const Eval& value)
1009 if constexpr (std::is_same_v<Value, Scalar>) {
1010 static_cast<void>(
this);
1011 return getValue(value);
1013 return this->extendEval(value);
1016 auto obtainN = [](
const auto& value)
1018 if constexpr (std::is_same_v<Value, Scalar>) {
1019 return getValue(value);
1024 const auto& fs = int_quants.fluidState();
1026 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
1027 const Value rs = obtain(fs.Rs());
1028 const Value rv = obtain(fs.Rv());
1031 std::vector<Value> b_perfcells(this->num_components_, 0.0);
1033 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1034 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1038 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1039 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
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));
1047 this->computePerfRate(pressure_cell,
1055 obtainN(this->segments_.density(seg)),
1064 template <
typename TypeTag>
1066 MultisegmentWell<TypeTag>::
1067 computeSegmentFluidProperties(
const Simulator& simulator, DeferredLogger& deferred_logger)
1076 EvalWell temperature;
1077 EvalWell saltConcentration;
1083 int pvt_region_index;
1086 const int cell_idx = this->well_cells_[0];
1087 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 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();
1094 this->segments_.computeFluidProperties(temperature,
1096 this->primary_variables_,
1101 template <
typename TypeTag>
1102 template<
class Value>
1104 MultisegmentWell<TypeTag>::
1105 getMobility(
const Simulator& simulator,
1107 std::vector<Value>& mob,
1108 DeferredLogger& deferred_logger)
const
1110 auto obtain = [
this](
const Eval& value)
1112 if constexpr (std::is_same_v<Value, Scalar>) {
1113 static_cast<void>(
this);
1114 return getValue(value);
1116 return this->extendEval(value);
1120 WellInterface<TypeTag>::getMobility(simulator, perf, mob, obtain, deferred_logger);
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());
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;
1142 template<
typename TypeTag>
1144 MultisegmentWell<TypeTag>::
1145 getRefDensity()
const
1147 return this->segments_.getRefDensity();
1150 template<
typename TypeTag>
1152 MultisegmentWell<TypeTag>::
1153 checkOperabilityUnderBHPLimit(
const WellState& ,
const Simulator& simulator, DeferredLogger& deferred_logger)
1155 const auto& summaryState = simulator.vanguard().summaryState();
1156 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1159 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1160 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1163 double total_ipr_mass_rate = 0.0;
1164 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1166 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1170 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1171 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1173 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1174 total_ipr_mass_rate += ipr_rate * rho;
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;
1181 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1185 std::vector<double> well_rates_bhp_limit;
1186 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1188 const double thp_limit = this->getTHPConstraint(summaryState);
1189 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1191 this->getRefDensity(),
1192 this->wellEcl().alq_value(summaryState),
1195 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1196 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1207 this->operability_status_.operable_under_only_bhp_limit =
true;
1208 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1214 template<
typename TypeTag>
1216 MultisegmentWell<TypeTag>::
1217 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
1222 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1223 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1225 const int nseg = this->numberOfSegments();
1226 std::vector<double> seg_dp(nseg, 0.0);
1227 for (
int seg = 0; seg < nseg; ++seg) {
1229 const double dp = this->getSegmentDp(seg,
1230 this->segments_.density(seg).value(),
1233 for (
const int perf : this->segments_.perforations()[seg]) {
1234 std::vector<Scalar> mob(this->num_components_, 0.0);
1237 getMobility(simulator, perf, mob, deferred_logger);
1239 const int cell_idx = this->well_cells_[perf];
1240 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
1241 const auto& fs = int_quantities.fluidState();
1243 const double perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1245 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1246 const double pressure_cell = this->getPerfCellPressure(fs).value();
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)) {
1254 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1255 b_perf[comp_idx] = fs.invB(phase).value();
1259 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1260 const double pressure_diff = pressure_cell - h_perf;
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());
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;
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();
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];
1290 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1291 ipr_a_perf[oil_comp_idx] += vap_oil_a;
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];
1296 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1297 ipr_b_perf[oil_comp_idx] += vap_oil_b;
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];
1308 template<
typename TypeTag>
1310 MultisegmentWell<TypeTag>::
1311 updateIPRImplicit(
const Simulator& simulator, WellState& well_state, DeferredLogger& deferred_logger)
1320 auto rates = well_state.well(this->index_of_well_).surface_rates;
1322 for (std::size_t p = 0; p < rates.size(); ++p) {
1323 zero_rates &= rates[p] == 0.0;
1325 auto& ws = well_state.well(this->index_of_well_);
1327 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
1328 deferred_logger.debug(msg);
1340 const auto& group_state = simulator.problem().wellModel().groupState();
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.);
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;
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);
1356 BVectorWell rhs(this->numberOfSegments());
1358 rhs[0][SPres] = -1.0;
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) {
1367 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
1369 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
1372 ws.production_cmode = cmode;
1375 template<
typename TypeTag>
1377 MultisegmentWell<TypeTag>::
1378 checkOperabilityUnderTHPLimit(
1379 const Simulator& simulator,
1380 const WellState& well_state,
1381 DeferredLogger& deferred_logger)
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);
1390 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1392 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1393 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
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);
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);
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 ");
1428 template<
typename TypeTag>
1430 MultisegmentWell<TypeTag>::
1431 iterateWellEqWithControl(
const Simulator& simulator,
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)
1439 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1441 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1445 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1450 std::vector<std::vector<Scalar> > residual_history;
1451 std::vector<double> measure_history;
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_) {
1462 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1464 BVectorWell dx_well;
1466 dx_well = this->linSys_.solve();
1468 catch(
const NumericalProblem& exp) {
1472 deferred_logger.problem(
"In MultisegmentWell::iterateWellEqWithControl for well "
1473 + this->name() +
": "+exp.what());
1477 if (it > this->param_.strict_inner_iter_wells_) {
1478 relax_convergence =
true;
1479 this->regularize_ =
true;
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()) {
1491 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
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_,
1504 bool is_oscillate =
false;
1505 bool is_stagnate =
false;
1507 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1511 if (is_oscillate || is_stagnate) {
1513 std::ostringstream sstr;
1514 if (relaxation_factor == min_relaxation_factor) {
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()) {
1522 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1523 deferred_logger.debug(sstr.str());
1530 const double reduction_mutliplier = 0.9;
1531 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1535 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1539 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1541 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1543 this->regularize_ =
true;
1544 deferred_logger.debug(sstr.str());
1546 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1547 initPrimaryVariablesEvaluation();
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());
1558 std::ostringstream sstr;
1559 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1560#define EXTRA_DEBUG_MSW 0
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) {
1569 sstr <<
" " << measure_history[i] <<
" \n";
1572#undef EXTRA_DEBUG_MSW
1573 deferred_logger.debug(sstr.str());
1580 template<
typename TypeTag>
1582 MultisegmentWell<TypeTag>::
1583 iterateWellEqWithSwitching(
const Simulator& simulator,
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 ,
1591 const bool fixed_status )
1593 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1597 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1602 std::vector<std::vector<Scalar> > residual_history;
1603 std::vector<double> measure_history;
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();
1617 const int min_its_after_switch = 3;
1618 int its_since_last_switch = min_its_after_switch;
1619 int switch_count= 0;
1621 const auto well_status_orig = this->wellStatus_;
1622 const auto operability_orig = this->operability_status_;
1623 auto well_status_cur = well_status_orig;
1625 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
1626 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
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;
1633 this->operability_status_.resetOperability();
1634 this->operability_status_.solvable =
true;
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);
1642 its_since_last_switch = 0;
1644 if (well_status_cur != this->wellStatus_) {
1645 well_status_cur = this->wellStatus_;
1648 if (!changed && final_check) {
1651 final_check =
false;
1655 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1657 const BVectorWell dx_well = this->linSys_.solve();
1659 if (it > this->param_.strict_inner_iter_wells_) {
1660 relax_convergence =
true;
1661 this->regularize_ =
true;
1664 const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1665 converged = report.converged();
1669 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1671 its_since_last_switch = min_its_after_switch;
1679 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1683 residual_history.push_back(residuals);
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_,
1693 bool is_oscillate =
false;
1694 bool is_stagnate =
false;
1696 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1700 if (is_oscillate || is_stagnate) {
1702 std::string message;
1703 if (relaxation_factor == min_relaxation_factor) {
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()) {
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);
1720 constexpr double reduction_mutliplier = 0.9;
1721 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1725 fmt::format_to(std::back_inserter(message),
" well {} observes stagnation in inner iteration {}\n", this->name(), it);
1728 fmt::format_to(std::back_inserter(message),
" well {} observes oscillation in inner iteration {}\n", this->name(), it);
1730 fmt::format_to(std::back_inserter(message),
" relaxation_factor is {} now\n", relaxation_factor);
1732 this->regularize_ =
true;
1733 deferred_logger.debug(message);
1736 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1737 initPrimaryVariablesEvaluation();
1741 if (allow_switching){
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;
1748 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
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_));
1757 deferred_logger.debug(message);
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);
1771 template<
typename TypeTag>
1773 MultisegmentWell<TypeTag>::
1774 assembleWellEqWithoutIteration(
const Simulator& simulator,
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)
1782 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1785 this->segments_.updateUpwindingSegments(this->primary_variables_);
1788 computeSegmentFluidProperties(simulator, deferred_logger);
1791 this->linSys_.clear();
1793 auto& ws = well_state.well(this->index_of_well_);
1794 ws.phase_mixing_rates.fill(0.0);
1801 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1803 const int nseg = this->numberOfSegments();
1805 for (
int seg = 0; seg < nseg; ++seg) {
1809 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg);
1814 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
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_);
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,
1831 this->well_efficiency_factor_;
1832 MultisegmentWellAssemble(*this).
1833 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
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,
1846 this->well_efficiency_factor_;
1847 MultisegmentWellAssemble(*this).
1848 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
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, 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);
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;
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();
1884 perf_press_state[perf] = perf_press.value();
1886 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1888 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1890 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1892 MultisegmentWellAssemble(*this).
1893 assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
1899 const auto& summaryState = simulator.vanguard().summaryState();
1900 const Schedule& schedule = simulator.vanguard().schedule();
1901 MultisegmentWellAssemble(*this).
1902 assembleControlEq(well_state,
1909 this->primary_variables_,
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);
1919 this->linSys_.createSolver();
1925 template<
typename TypeTag>
1927 MultisegmentWell<TypeTag>::
1928 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1930 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1934 template<
typename TypeTag>
1936 MultisegmentWell<TypeTag>::
1937 allDrawDownWrongDirection(
const Simulator& simulator)
const
1939 bool all_drawdown_wrong_direction =
true;
1940 const int nseg = this->numberOfSegments();
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]) {
1946 const int cell_idx = this->well_cells_[perf];
1947 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1948 const auto& fs = intQuants.fluidState();
1951 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1953 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1955 const double pressure_cell = this->getPerfCellPressure(fs).value();
1956 const double perf_press = pressure_cell - cell_perf_press_diff;
1959 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1964 if ( (drawdown < 0. && this->isInjector()) ||
1965 (drawdown > 0. && this->isProducer()) ) {
1966 all_drawdown_wrong_direction =
false;
1972 return all_drawdown_wrong_direction;
1978 template<
typename TypeTag>
1980 MultisegmentWell<TypeTag>::
1981 updateWaterThroughput(
const double , WellState& )
const
1989 template<
typename TypeTag>
1990 typename MultisegmentWell<TypeTag>::EvalWell
1991 MultisegmentWell<TypeTag>::
1992 getSegmentSurfaceVolume(
const Simulator& simulator,
const int seg_idx)
const
1994 EvalWell temperature;
1995 EvalWell saltConcentration;
1996 int pvt_region_index;
2000 const int cell_idx = this->well_cells_[0];
2001 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 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();
2008 return this->segments_.getSurfaceVolume(temperature,
2010 this->primary_variables_,
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
2024 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
2027 this->getALQ(well_state),
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
2042 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2048 std::vector<double> rates(3);
2049 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2053 auto bhpAtLimit = WellBhpThpCalculator(*this).
2054 computeBhpAtThpLimitProd(frates,
2056 this->maxPerfPress(simulator),
2057 this->getRefDensity(),
2059 this->getTHPConstraint(summary_state),
2065 auto fratesIter = [
this, &simulator, &deferred_logger](
const double bhp) {
2069 std::vector<double> rates(3);
2070 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2074 return WellBhpThpCalculator(*this).
2075 computeBhpAtThpLimitProd(fratesIter,
2077 this->maxPerfPress(simulator),
2078 this->getRefDensity(),
2080 this->getTHPConstraint(summary_state),
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
2092 auto frates = [
this, &simulator, &deferred_logger](
const double bhp) {
2098 std::vector<double> rates(3);
2099 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2103 auto bhpAtLimit = WellBhpThpCalculator(*this).
2104 computeBhpAtThpLimitInj(frates,
2106 this->getRefDensity(),
2115 auto fratesIter = [
this, &simulator, &deferred_logger](
const double bhp) {
2119 std::vector<double> rates(3);
2120 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2124 return WellBhpThpCalculator(*this).
2125 computeBhpAtThpLimitInj(fratesIter,
2127 this->getRefDensity(),
2138 template<
typename TypeTag>
2140 MultisegmentWell<TypeTag>::
2141 maxPerfPress(
const Simulator& simulator)
const
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, 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);
2154 return max_pressure;
2161 template<
typename TypeTag>
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) {
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, 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;
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];
2196 template <
typename TypeTag>
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);
2214 template <
typename TypeTag>
2216 MultisegmentWell<TypeTag>::
2217 setPrimaryVars(std::vector<double>::const_iterator it)
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);
2227 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 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