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,Scalar>&>(*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>
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,
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_,
217 template <
typename TypeTag>
220 apply(
const BVector& x, BVector&
Ax)
const
222 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
226 if (this->param_.matrix_add_well_contributions_) {
238 template <
typename TypeTag>
243 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
252 template <
typename TypeTag>
260 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
265 this->linSys_.recoverSolutionWell(x,
xw);
273 template <
typename TypeTag>
278 std::vector<double>& well_potentials,
282 this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
288 debug_cost_counter_ = 0;
290 const auto& summaryState = ebosSimulator.vanguard().summaryState();
292 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials,
deferred_logger);
294 well_potentials = computeWellPotentialWithTHP(
297 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
298 + this->name() +
": " + std::to_string(debug_cost_counter_));
300 this->checkNegativeWellPotentials(well_potentials,
301 this->param_.check_well_operability_,
308 template<
typename TypeTag>
315 if (this->well_ecl_.isInjector()) {
316 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
319 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
320 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
324 template<
typename TypeTag>
326 MultisegmentWell<TypeTag>::
327 computeWellRatesWithBhp(
const Simulator& ebosSimulator,
329 std::vector<double>& well_flux,
330 DeferredLogger& deferred_logger)
const
333 const int np = this->number_of_phases_;
335 well_flux.resize(np, 0.0);
336 const bool allow_cf = this->getAllowCrossFlow();
337 const int nseg = this->numberOfSegments();
338 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
339 const auto& ws = well_state.well(this->indexOfWell());
340 auto segments_copy = ws.segments;
341 segments_copy.scale_pressure(bhp);
342 const auto& segment_pressure = segments_copy.pressure;
343 for (
int seg = 0; seg < nseg; ++seg) {
344 for (
const int perf : this->segments_.perforations()[seg]) {
345 const int cell_idx = this->well_cells_[perf];
346 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
348 std::vector<Scalar> mob(this->num_components_, 0.);
349 getMobility(ebosSimulator, perf, mob, deferred_logger);
350 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
351 const double Tw = this->well_index_[perf] * trans_mult;
353 const Scalar seg_pressure = segment_pressure[seg];
354 std::vector<Scalar> cq_s(this->num_components_, 0.);
355 Scalar perf_press = 0.0;
356 PerforationRates perf_rates;
357 computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
358 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
360 for(
int p = 0; p < np; ++p) {
361 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
365 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
369 template<
typename TypeTag>
371 MultisegmentWell<TypeTag>::
372 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
374 std::vector<double>& well_flux,
375 DeferredLogger& deferred_logger)
const
379 MultisegmentWell<TypeTag> well_copy(*
this);
380 well_copy.debug_cost_counter_ = 0;
383 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
384 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
385 auto& ws = well_state_copy.well(this->index_of_well_);
388 const auto& summary_state = ebosSimulator.vanguard().summaryState();
389 auto inj_controls = well_copy.well_ecl_.isInjector()
390 ? well_copy.well_ecl_.injectionControls(summary_state)
391 : Well::InjectionControls(0);
392 auto prod_controls = well_copy.well_ecl_.isProducer()
393 ? well_copy.well_ecl_.productionControls(summary_state) :
394 Well::ProductionControls(0);
397 if (well_copy.well_ecl_.isInjector()) {
398 inj_controls.bhp_limit = bhp;
399 ws.injection_cmode = Well::InjectorCMode::BHP;
401 prod_controls.bhp_limit = bhp;
402 ws.production_cmode = Well::ProducerCMode::BHP;
405 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
408 const int np = this->number_of_phases_;
410 for (
int phase = 0; phase < np; ++phase){
411 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
414 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
415 for (
int phase = 0; phase < np; ++phase) {
416 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
419 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
420 this->segments_.perforations(),
423 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
424 const double dt = ebosSimulator.timeStepSize();
426 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
431 well_flux.resize(np, 0.0);
432 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
433 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
434 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
436 debug_cost_counter_ += well_copy.debug_cost_counter_;
441 template<
typename TypeTag>
443 MultisegmentWell<TypeTag>::
444 computeWellPotentialWithTHP(
445 const WellState& well_state,
446 const Simulator& ebos_simulator,
447 DeferredLogger& deferred_logger)
const
449 std::vector<double> potentials(this->number_of_phases_, 0.0);
450 const auto& summary_state = ebos_simulator.vanguard().summaryState();
452 const auto& well = this->well_ecl_;
453 if (well.isInjector()){
454 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
455 if (bhp_at_thp_limit) {
456 const auto& controls = well.injectionControls(summary_state);
457 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
458 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
459 deferred_logger.debug(
"Converged thp based potential calculation for well "
460 + this->name() +
", at bhp = " + std::to_string(bhp));
462 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
463 "Failed in getting converged thp based potential calculation for well "
464 + this->name() +
". Instead the bhp based value is used");
465 const auto& controls = well.injectionControls(summary_state);
466 const double bhp = controls.bhp_limit;
467 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
470 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
471 well_state, ebos_simulator, summary_state, deferred_logger);
472 if (bhp_at_thp_limit) {
473 const auto& controls = well.productionControls(summary_state);
474 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
475 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
476 deferred_logger.debug(
"Converged thp based potential calculation for well "
477 + this->name() +
", at bhp = " + std::to_string(bhp));
479 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
480 "Failed in getting converged thp based potential calculation for well "
481 + this->name() +
". Instead the bhp based value is used");
482 const auto& controls = well.productionControls(summary_state);
483 const double bhp = controls.bhp_limit;
484 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
493 template <
typename TypeTag>
495 MultisegmentWell<TypeTag>::
496 solveEqAndUpdateWellState(
const SummaryState& summary_state,
497 WellState& well_state,
498 DeferredLogger& deferred_logger)
500 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
504 const BVectorWell dx_well = this->linSys_.solve();
506 updateWellState(summary_state, dx_well, well_state, deferred_logger);
513 template <
typename TypeTag>
515 MultisegmentWell<TypeTag>::
516 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
518 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
520 std::vector<double> kr(this->number_of_phases_, 0.0);
521 std::vector<double> density(this->number_of_phases_, 0.0);
523 const int cell_idx = this->well_cells_[perf];
524 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
525 const auto& fs = intQuants.fluidState();
530 if (pu.phase_used[Water]) {
531 const int water_pos = pu.phase_pos[Water];
532 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
533 sum_kr += kr[water_pos];
534 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
537 if (pu.phase_used[Oil]) {
538 const int oil_pos = pu.phase_pos[Oil];
539 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
540 sum_kr += kr[oil_pos];
541 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
544 if (pu.phase_used[Gas]) {
545 const int gas_pos = pu.phase_pos[Gas];
546 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
547 sum_kr += kr[gas_pos];
548 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
551 assert(sum_kr != 0.);
554 double average_density = 0.;
555 for (
int p = 0; p < this->number_of_phases_; ++p) {
556 average_density += kr[p] * density[p];
558 average_density /= sum_kr;
560 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
568 template <
typename TypeTag>
570 MultisegmentWell<TypeTag>::
571 computeInitialSegmentFluids(
const Simulator& ebos_simulator)
573 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
575 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
576 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
577 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
586 template <
typename TypeTag>
588 MultisegmentWell<TypeTag>::
589 updateWellState(
const SummaryState& summary_state,
590 const BVectorWell& dwells,
591 WellState& well_state,
592 DeferredLogger& deferred_logger,
593 const double relaxation_factor)
595 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
597 const double dFLimit = this->param_.dwell_fraction_max_;
598 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
599 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
600 this->primary_variables_.updateNewton(dwells,
603 stop_or_zero_rate_target,
604 max_pressure_change);
606 this->primary_variables_.copyToWellState(*
this, getRefDensity(), stop_or_zero_rate_target,
607 well_state, summary_state, deferred_logger);
610 auto& ws = well_state.well(this->index_of_well_);
611 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
614 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
621 template <
typename TypeTag>
623 MultisegmentWell<TypeTag>::
624 calculateExplicitQuantities(
const Simulator& ebosSimulator,
625 const WellState& well_state,
626 DeferredLogger& deferred_logger)
628 const auto& summary_state = ebosSimulator.vanguard().summaryState();
629 updatePrimaryVariables(summary_state, well_state, deferred_logger);
630 initPrimaryVariablesEvaluation();
631 computePerfCellPressDiffs(ebosSimulator);
632 computeInitialSegmentFluids(ebosSimulator);
639 template<
typename TypeTag>
641 MultisegmentWell<TypeTag>::
642 updateProductivityIndex(
const Simulator& ebosSimulator,
643 const WellProdIndexCalculator& wellPICalc,
644 WellState& well_state,
645 DeferredLogger& deferred_logger)
const
647 auto fluidState = [&ebosSimulator,
this](
const int perf)
649 const auto cell_idx = this->well_cells_[perf];
650 return ebosSimulator.model()
651 .intensiveQuantities(cell_idx, 0).fluidState();
654 const int np = this->number_of_phases_;
655 auto setToZero = [np](
double* x) ->
void
657 std::fill_n(x, np, 0.0);
660 auto addVector = [np](
const double* src,
double* dest) ->
void
662 std::transform(src, src + np, dest, dest, std::plus<>{});
665 auto& ws = well_state.well(this->index_of_well_);
666 auto& perf_data = ws.perf_data;
667 auto* connPI = perf_data.prod_index.data();
668 auto* wellPI = ws.productivity_index.data();
672 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
673 auto subsetPerfID = 0;
675 for (
const auto& perf : *this->perf_data_){
676 auto allPerfID = perf.ecl_index;
678 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
680 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
683 std::vector<Scalar> mob(this->num_components_, 0.0);
684 getMobility(ebosSimulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
686 const auto& fs = fluidState(subsetPerfID);
689 if (this->isInjector()) {
690 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
691 mob, connPI, deferred_logger);
694 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
697 addVector(connPI, wellPI);
703 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
704 "Internal logic error in processing connections for PI/II");
711 template<
typename TypeTag>
713 MultisegmentWell<TypeTag>::
714 connectionDensity(
const int globalConnIdx,
715 [[maybe_unused]]
const int openConnIdx)
const
720 const auto segNum = this->wellEcl()
721 .getConnections()[globalConnIdx].segment();
723 const auto segIdx = this->wellEcl()
724 .getSegments().segmentNumberToIndex(segNum);
726 return this->segments_.density(segIdx).value();
733 template<
typename TypeTag>
735 MultisegmentWell<TypeTag>::
736 addWellContributions(SparseMatrixAdapter& jacobian)
const
738 this->linSys_.extract(jacobian);
742 template<
typename TypeTag>
744 MultisegmentWell<TypeTag>::
745 addWellPressureEquations(PressureMatrix& jacobian,
746 const BVector& weights,
747 const int pressureVarIndex,
748 const bool use_well_weights,
749 const WellState& well_state)
const
752 this->linSys_.extractCPRPressureMatrix(jacobian,
762 template<
typename TypeTag>
763 template<
class Value>
765 MultisegmentWell<TypeTag>::
766 computePerfRate(
const Value& pressure_cell,
769 const std::vector<Value>& b_perfcells,
770 const std::vector<Value>& mob_perfcells,
773 const Value& segment_pressure,
774 const Value& segment_density,
775 const bool& allow_cf,
776 const std::vector<Value>& cmix_s,
777 std::vector<Value>& cq_s,
779 PerforationRates& perf_rates,
780 DeferredLogger& deferred_logger)
const
783 const Value perf_seg_press_diff = this->gravity() * segment_density *
784 this->segments_.perforation_depth_diff(perf);
786 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
790 perf_press = segment_pressure + perf_seg_press_diff;
793 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
796 const Value drawdown = cell_press_at_perf - perf_press;
799 if (drawdown > 0.0) {
801 if (!allow_cf && this->isInjector()) {
806 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
807 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
808 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
811 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
812 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
813 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
814 const Value cq_s_oil = cq_s[oilCompIdx];
815 const Value cq_s_gas = cq_s[gasCompIdx];
816 cq_s[gasCompIdx] += rs * cq_s_oil;
817 cq_s[oilCompIdx] += rv * cq_s_gas;
821 if (!allow_cf && this->isProducer()) {
826 Value total_mob = mob_perfcells[0];
827 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
828 total_mob += mob_perfcells[comp_idx];
832 const Value cqt_i = - Tw * (total_mob * drawdown);
835 Value volume_ratio = 0.0;
836 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
837 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
838 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
841 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
842 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
843 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
848 const Value d = 1.0 - rv * rs;
850 if (getValue(d) == 0.0) {
851 OPM_DEFLOG_THROW(NumericalProblem,
852 fmt::format(
"Zero d value obtained for well {} "
853 "during flux calculation with rs {} and rv {}",
854 this->name(), rs, rv),
858 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
859 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
861 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
862 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
864 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
865 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
866 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
868 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
869 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
870 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
874 Value cqt_is = cqt_i / volume_ratio;
875 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
876 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
881 if (this->isProducer()) {
882 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
883 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
884 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
893 const double d = 1.0 - getValue(rv) * getValue(rs);
896 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
899 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
904 template <
typename TypeTag>
905 template<
class Value>
907 MultisegmentWell<TypeTag>::
908 computePerfRate(
const IntensiveQuantities& int_quants,
909 const std::vector<Value>& mob_perfcells,
913 const Value& segment_pressure,
914 const bool& allow_cf,
915 std::vector<Value>& cq_s,
917 PerforationRates& perf_rates,
918 DeferredLogger& deferred_logger)
const
921 auto obtain = [
this](
const Eval& value)
923 if constexpr (std::is_same_v<Value, Scalar>) {
924 static_cast<void>(
this);
925 return getValue(value);
927 return this->extendEval(value);
930 auto obtainN = [](
const auto& value)
932 if constexpr (std::is_same_v<Value, Scalar>) {
933 return getValue(value);
938 const auto& fs = int_quants.fluidState();
940 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
941 const Value rs = obtain(fs.Rs());
942 const Value rv = obtain(fs.Rv());
945 std::vector<Value> b_perfcells(this->num_components_, 0.0);
947 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
948 if (!FluidSystem::phaseIsActive(phaseIdx)) {
952 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
953 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
956 std::vector<Value> cmix_s(this->numComponents(), 0.0);
957 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
958 cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
961 this->computePerfRate(pressure_cell,
969 obtainN(this->segments_.density(seg)),
978 template <
typename TypeTag>
980 MultisegmentWell<TypeTag>::
981 computeSegmentFluidProperties(
const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
990 EvalWell temperature;
991 EvalWell saltConcentration;
997 int pvt_region_index;
1000 const int cell_idx = this->well_cells_[0];
1001 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1002 const auto& fs = intQuants.fluidState();
1003 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1004 saltConcentration = this->extendEval(fs.saltConcentration());
1005 pvt_region_index = fs.pvtRegionIndex();
1008 this->segments_.computeFluidProperties(temperature,
1010 this->primary_variables_,
1015 template <
typename TypeTag>
1016 template<
class Value>
1018 MultisegmentWell<TypeTag>::
1019 getMobility(
const Simulator& ebosSimulator,
1021 std::vector<Value>& mob,
1022 DeferredLogger& deferred_logger)
const
1024 auto obtain = [
this](
const Eval& value)
1026 if constexpr (std::is_same_v<Value, Scalar>) {
1027 static_cast<void>(
this);
1028 return getValue(value);
1030 return this->extendEval(value);
1034 WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
1036 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
1037 const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
1038 const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
1039 const int seg = this->segmentNumberToIndex(con.segment());
1043 const double segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
1044 const double perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
1045 * this->segments_.perforation_depth_diff(perf);
1046 const double perf_press = segment_pres + perf_seg_press_diff;
1047 const double multiplier = this->getInjMult(perf, segment_pres, perf_press);
1048 for (std::size_t i = 0; i < mob.size(); ++i) {
1049 mob[i] *= multiplier;
1056 template<
typename TypeTag>
1058 MultisegmentWell<TypeTag>::
1059 getRefDensity()
const
1061 return this->segments_.getRefDensity();
1064 template<
typename TypeTag>
1066 MultisegmentWell<TypeTag>::
1067 checkOperabilityUnderBHPLimit(
const WellState& ,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1069 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1070 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1073 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1074 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1077 double total_ipr_mass_rate = 0.0;
1078 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1080 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1084 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1085 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1087 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1088 total_ipr_mass_rate += ipr_rate * rho;
1090 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1091 this->operability_status_.operable_under_only_bhp_limit =
false;
1095 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1099 std::vector<double> well_rates_bhp_limit;
1100 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1102 const double thp_limit = this->getTHPConstraint(summaryState);
1103 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1105 this->getRefDensity(),
1106 this->wellEcl().alq_value(),
1109 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1110 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1121 this->operability_status_.operable_under_only_bhp_limit =
true;
1122 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1128 template<
typename TypeTag>
1130 MultisegmentWell<TypeTag>::
1131 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1136 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1137 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1139 const int nseg = this->numberOfSegments();
1140 std::vector<double> seg_dp(nseg, 0.0);
1141 for (
int seg = 0; seg < nseg; ++seg) {
1143 const double dp = this->getSegmentDp(seg,
1144 this->segments_.density(seg).value(),
1147 for (
const int perf : this->segments_.perforations()[seg]) {
1148 std::vector<Scalar> mob(this->num_components_, 0.0);
1151 getMobility(ebos_simulator, perf, mob, deferred_logger);
1153 const int cell_idx = this->well_cells_[perf];
1154 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1155 const auto& fs = int_quantities.fluidState();
1157 const double perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1159 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1160 const double pressure_cell = this->getPerfCellPressure(fs).value();
1163 std::vector<double> b_perf(this->num_components_);
1164 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1165 if (!FluidSystem::phaseIsActive(phase)) {
1168 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1169 b_perf[comp_idx] = fs.invB(phase).value();
1173 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1174 const double pressure_diff = pressure_cell - h_perf;
1177 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1178 deferred_logger.debug(
"CROSSFLOW_IPR",
1179 "cross flow found when updateIPR for well " + this->name());
1183 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1185 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1186 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1187 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1188 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1189 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1190 ipr_b_perf[comp_idx] += tw_mob;
1194 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1195 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1196 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1197 const double rs = (fs.Rs()).value();
1198 const double rv = (fs.Rv()).value();
1200 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1201 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1203 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1204 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1206 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1207 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1209 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1210 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1213 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1214 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1215 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1221 template<
typename TypeTag>
1223 MultisegmentWell<TypeTag>::
1224 checkOperabilityUnderTHPLimit(
1225 const Simulator& ebos_simulator,
1226 const WellState& well_state,
1227 DeferredLogger& deferred_logger)
1229 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1230 const auto obtain_bhp = this->isProducer()
1231 ? computeBhpAtThpLimitProd(
1232 well_state, ebos_simulator, summaryState, deferred_logger)
1233 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1236 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1238 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1239 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1241 const double thp_limit = this->getTHPConstraint(summaryState);
1242 if (this->isProducer() && *obtain_bhp < thp_limit) {
1243 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1244 +
" bars is SMALLER than thp limit "
1245 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1246 +
" bars as a producer for well " + this->name();
1247 deferred_logger.debug(msg);
1249 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1250 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1251 +
" bars is LARGER than thp limit "
1252 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1253 +
" bars as a injector for well " + this->name();
1254 deferred_logger.debug(msg);
1259 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1260 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1261 if (!this->wellIsStopped()) {
1262 const double thp_limit = this->getTHPConstraint(summaryState);
1263 deferred_logger.debug(
" could not find bhp value at thp limit "
1264 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1265 +
" bar for well " + this->name() +
", the well might need to be closed ");
1274 template<
typename TypeTag>
1276 MultisegmentWell<TypeTag>::
1277 iterateWellEqWithControl(
const Simulator& ebosSimulator,
1279 const Well::InjectionControls& inj_controls,
1280 const Well::ProductionControls& prod_controls,
1281 WellState& well_state,
1282 const GroupState& group_state,
1283 DeferredLogger& deferred_logger)
1285 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1287 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1291 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1296 std::vector<std::vector<Scalar> > residual_history;
1297 std::vector<double> measure_history;
1300 double relaxation_factor = 1.;
1301 const double min_relaxation_factor = 0.6;
1302 bool converged =
false;
1303 int stagnate_count = 0;
1304 bool relax_convergence =
false;
1305 this->regularize_ =
false;
1306 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1308 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1310 const BVectorWell dx_well = this->linSys_.solve();
1312 if (it > this->param_.strict_inner_iter_wells_) {
1313 relax_convergence =
true;
1314 this->regularize_ =
true;
1317 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1318 const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1319 if (report.converged()) {
1326 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1330 residual_history.push_back(residuals);
1331 measure_history.push_back(this->getResidualMeasureValue(well_state,
1332 residual_history[it],
1333 this->param_.tolerance_wells_,
1334 this->param_.tolerance_pressure_ms_wells_,
1339 bool is_oscillate =
false;
1340 bool is_stagnate =
false;
1342 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1346 if (is_oscillate || is_stagnate) {
1348 std::ostringstream sstr;
1349 if (relaxation_factor == min_relaxation_factor) {
1352 if (stagnate_count == 6) {
1353 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1354 const auto reportStag = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger,
true);
1355 if (reportStag.converged()) {
1357 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1358 deferred_logger.debug(sstr.str());
1365 const double reduction_mutliplier = 0.9;
1366 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1370 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1374 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1376 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1378 this->regularize_ =
true;
1379 deferred_logger.debug(sstr.str());
1381 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1382 initPrimaryVariablesEvaluation();
1387 std::ostringstream sstr;
1388 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1389 if (relax_convergence)
1390 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1391 deferred_logger.debug(sstr.str());
1393 std::ostringstream sstr;
1394 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1395#define EXTRA_DEBUG_MSW 0
1397 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1398 for (
int i = 0; i < it; ++i) {
1399 const auto& residual = residual_history[i];
1400 sstr <<
" residual at " << i <<
"th iteration ";
1401 for (
const auto& res : residual) {
1404 sstr <<
" " << measure_history[i] <<
" \n";
1407#undef EXTRA_DEBUG_MSW
1408 deferred_logger.debug(sstr.str());
1415 template<
typename TypeTag>
1417 MultisegmentWell<TypeTag>::
1418 iterateWellEqWithSwitching(
const Simulator& ebosSimulator,
1420 const Well::InjectionControls& inj_controls,
1421 const Well::ProductionControls& prod_controls,
1422 WellState& well_state,
1423 const GroupState& group_state,
1424 DeferredLogger& deferred_logger)
1428 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1432 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1437 std::vector<std::vector<Scalar> > residual_history;
1438 std::vector<double> measure_history;
1441 double relaxation_factor = 1.;
1442 const double min_relaxation_factor = 0.6;
1443 bool converged =
false;
1444 [[maybe_unused]]
int stagnate_count = 0;
1445 bool relax_convergence =
false;
1446 this->regularize_ =
false;
1449 const int min_its_after_switch = 2;
1450 int its_since_last_switch = min_its_after_switch;
1451 int switch_count= 0;
1452 const auto well_status = this->wellStatus_;
1453 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1454 const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
1455 bool changed =
false;
1456 bool final_check =
false;
1458 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1459 its_since_last_switch++;
1460 if (its_since_last_switch >= min_its_after_switch){
1461 const double wqTotal = this->primary_variables_.getWQTotal().value();
1462 changed = this->updateWellControlAndStatusLocalIteration (ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger);
1464 its_since_last_switch = 0;
1467 if (!changed && final_check) {
1470 final_check =
false;
1474 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1476 const BVectorWell dx_well = this->linSys_.solve();
1478 if (it > this->param_.strict_inner_iter_wells_) {
1479 relax_convergence =
true;
1480 this->regularize_ =
true;
1483 const auto report = getWellConvergence(summary_state, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1484 converged = report.converged();
1488 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1490 its_since_last_switch = min_its_after_switch;
1498 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1502 residual_history.push_back(residuals);
1506 measure_history.push_back(this->getResidualMeasureValue(well_state,
1507 residual_history[it],
1508 this->param_.tolerance_wells_,
1509 this->param_.tolerance_pressure_ms_wells_,
1512 bool is_oscillate =
false;
1513 bool is_stagnate =
false;
1515 this->detectOscillations(measure_history, is_oscillate, is_stagnate);
1519 if (is_oscillate || is_stagnate) {
1521 std::string message;
1522 if (relaxation_factor == min_relaxation_factor) {
1525 fmt::format_to(std::back_inserter(message),
" Well {} observes severe stagnation and/or oscillation."
1526 " We relax the tolerance and check for convergence. \n", this->name());
1527 const auto reportStag = getWellConvergence(summary_state, well_state, Base::B_avg_,
1528 deferred_logger,
true);
1529 if (reportStag.converged()) {
1531 fmt::format_to(std::back_inserter(message),
" Well {} manages to get converged with relaxed tolerances in {} inner iterations", this->name(), it);
1532 deferred_logger.debug(message);
1539 constexpr double reduction_mutliplier = 0.9;
1540 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1544 fmt::format_to(std::back_inserter(message),
" well {} observes stagnation in inner iteration {}\n", this->name(), it);
1547 fmt::format_to(std::back_inserter(message),
" well {} observes oscillation in inner iteration {}\n", this->name(), it);
1549 fmt::format_to(std::back_inserter(message),
" relaxation_factor is {} now\n", relaxation_factor);
1551 this->regularize_ =
true;
1552 deferred_logger.debug(message);
1555 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1556 initPrimaryVariablesEvaluation();
1560 if (allow_switching){
1562 const bool is_stopped = this->wellIsStopped();
1563 if (this->wellHasTHPConstraints(summary_state)){
1564 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
1565 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
1567 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
1571 this->wellStatus_ = well_status;
1573 std::string message = fmt::format(
" Well {} converged in {} inner iterations ("
1574 "{} control/status switches).", this->name(), it, switch_count);
1575 if (relax_convergence) {
1576 message.append(fmt::format(
" (A relaxed tolerance was used after {} iterations)",
1577 this->param_.strict_inner_iter_wells_));
1579 deferred_logger.debug(message);
1581 const std::string message = fmt::format(
" Well {} did not converged in {} inner iterations ("
1582 "{} control/status switches).", this->name(), it, switch_count);
1583 deferred_logger.debug(message);
1590 template<
typename TypeTag>
1592 MultisegmentWell<TypeTag>::
1593 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
1595 const Well::InjectionControls& inj_controls,
1596 const Well::ProductionControls& prod_controls,
1597 WellState& well_state,
1598 const GroupState& group_state,
1599 DeferredLogger& deferred_logger)
1601 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1604 this->segments_.updateUpwindingSegments(this->primary_variables_);
1607 computeSegmentFluidProperties(ebosSimulator, deferred_logger);
1610 this->linSys_.clear();
1612 auto& ws = well_state.well(this->index_of_well_);
1613 ws.phase_mixing_rates.fill(0.0);
1620 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1622 const int nseg = this->numberOfSegments();
1624 for (
int seg = 0; seg < nseg; ++seg) {
1628 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1633 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1635 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1636 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1637 - segment_fluid_initial_[seg][comp_idx]) / dt;
1638 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1639 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1644 const int seg_upwind = this->segments_.upwinding_segment(seg);
1645 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1646 const EvalWell segment_rate =
1647 this->primary_variables_.getSegmentRateUpwinding(seg,
1650 this->well_efficiency_factor_;
1651 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1652 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1658 for (
const int inlet : this->segments_.inlets()[seg]) {
1659 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1660 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1661 const EvalWell inlet_rate =
1662 this->primary_variables_.getSegmentRateUpwinding(inlet,
1665 this->well_efficiency_factor_;
1666 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1667 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1673 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1674 auto& perf_data = ws.perf_data;
1675 auto& perf_rates = perf_data.phase_rates;
1676 auto& perf_press_state = perf_data.pressure;
1677 for (
const int perf : this->segments_.perforations()[seg]) {
1678 const int cell_idx = this->well_cells_[perf];
1679 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1680 std::vector<EvalWell> mob(this->num_components_, 0.0);
1681 getMobility(ebosSimulator, perf, mob, deferred_logger);
1682 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1683 const double Tw = this->well_index_[perf] * trans_mult;
1684 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1685 EvalWell perf_press;
1686 PerforationRates perfRates;
1687 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
1688 allow_cf, cq_s, perf_press, perfRates, deferred_logger);
1691 if (this->isProducer()) {
1692 ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas;
1693 ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil;
1697 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1698 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1700 perf_press_state[perf] = perf_press.value();
1702 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1704 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1706 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1708 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1709 assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
1715 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1716 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1717 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1718 assembleControlEq(well_state,
1725 this->primary_variables_,
1729 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1730 this->assemblePressureEq(seg, unit_system, well_state, this->param_.use_average_density_ms_wells_, deferred_logger);
1734 this->linSys_.createSolver();
1740 template<
typename TypeTag>
1742 MultisegmentWell<TypeTag>::
1743 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1745 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1749 template<
typename TypeTag>
1751 MultisegmentWell<TypeTag>::
1752 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1754 bool all_drawdown_wrong_direction =
true;
1755 const int nseg = this->numberOfSegments();
1757 for (
int seg = 0; seg < nseg; ++seg) {
1758 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1759 for (
const int perf : this->segments_.perforations()[seg]) {
1761 const int cell_idx = this->well_cells_[perf];
1762 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1763 const auto& fs = intQuants.fluidState();
1766 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1768 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1770 const double pressure_cell = this->getPerfCellPressure(fs).value();
1771 const double perf_press = pressure_cell - cell_perf_press_diff;
1774 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1779 if ( (drawdown < 0. && this->isInjector()) ||
1780 (drawdown > 0. && this->isProducer()) ) {
1781 all_drawdown_wrong_direction =
false;
1787 return all_drawdown_wrong_direction;
1793 template<
typename TypeTag>
1795 MultisegmentWell<TypeTag>::
1796 updateWaterThroughput(
const double , WellState& )
const
1804 template<
typename TypeTag>
1805 typename MultisegmentWell<TypeTag>::EvalWell
1806 MultisegmentWell<TypeTag>::
1807 getSegmentSurfaceVolume(
const Simulator& ebos_simulator,
const int seg_idx)
const
1809 EvalWell temperature;
1810 EvalWell saltConcentration;
1811 int pvt_region_index;
1815 const int cell_idx = this->well_cells_[0];
1816 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1817 const auto& fs = intQuants.fluidState();
1818 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1819 saltConcentration = this->extendEval(fs.saltConcentration());
1820 pvt_region_index = fs.pvtRegionIndex();
1823 return this->segments_.getSurfaceVolume(temperature,
1825 this->primary_variables_,
1831 template<
typename TypeTag>
1832 std::optional<double>
1833 MultisegmentWell<TypeTag>::
1834 computeBhpAtThpLimitProd(
const WellState& well_state,
1835 const Simulator& ebos_simulator,
1836 const SummaryState& summary_state,
1837 DeferredLogger& deferred_logger)
const
1839 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
1842 this->getALQ(well_state),
1848 template<
typename TypeTag>
1849 std::optional<double>
1850 MultisegmentWell<TypeTag>::
1851 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
1852 const SummaryState& summary_state,
1853 const double alq_value,
1854 DeferredLogger& deferred_logger)
const
1857 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1863 std::vector<double> rates(3);
1864 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1868 auto bhpAtLimit = WellBhpThpCalculator(*this).
1869 computeBhpAtThpLimitProd(frates,
1871 this->maxPerfPress(ebos_simulator),
1872 this->getRefDensity(),
1874 this->getTHPConstraint(summary_state),
1880 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1884 std::vector<double> rates(3);
1885 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1889 return WellBhpThpCalculator(*this).
1890 computeBhpAtThpLimitProd(fratesIter,
1892 this->maxPerfPress(ebos_simulator),
1893 this->getRefDensity(),
1895 this->getTHPConstraint(summary_state),
1899 template<
typename TypeTag>
1900 std::optional<double>
1901 MultisegmentWell<TypeTag>::
1902 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
1903 const SummaryState& summary_state,
1904 DeferredLogger& deferred_logger)
const
1907 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1913 std::vector<double> rates(3);
1914 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1918 auto bhpAtLimit = WellBhpThpCalculator(*this).
1919 computeBhpAtThpLimitInj(frates,
1921 this->getRefDensity(),
1930 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1934 std::vector<double> rates(3);
1935 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1939 return WellBhpThpCalculator(*this).
1940 computeBhpAtThpLimitInj(fratesIter,
1942 this->getRefDensity(),
1953 template<
typename TypeTag>
1955 MultisegmentWell<TypeTag>::
1956 maxPerfPress(
const Simulator& ebos_simulator)
const
1958 double max_pressure = 0.0;
1959 const int nseg = this->numberOfSegments();
1960 for (
int seg = 0; seg < nseg; ++seg) {
1961 for (
const int perf : this->segments_.perforations()[seg]) {
1962 const int cell_idx = this->well_cells_[perf];
1963 const auto& int_quants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1964 const auto& fs = int_quants.fluidState();
1965 double pressure_cell = this->getPerfCellPressure(fs).value();
1966 max_pressure = std::max(max_pressure, pressure_cell);
1969 return max_pressure;
1976 template<
typename TypeTag>
1983 std::vector<Scalar>
well_q_s(this->num_components_, 0.0);
1984 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1985 const int nseg = this->numberOfSegments();
1989 for (
const int perf : this->segments_.perforations()[
seg]) {
1992 std::vector<Scalar>
mob(this->num_components_, 0.0);
1996 std::vector<Scalar>
cq_s(this->num_components_, 0.0);
2001 for (
int comp = 0;
comp < this->num_components_; ++
comp) {
2010 template <
typename TypeTag>
2015 const int num_seg = this->numberOfSegments();
2016 constexpr int num_eq = MSWEval::numWellEq;
2019 const auto& pv = this->primary_variables_.value(
ii);
2028 template <
typename TypeTag>
2030 MultisegmentWell<TypeTag>::
2031 setPrimaryVars(std::vector<double>::const_iterator it)
2033 const int num_seg = this->numberOfSegments();
2034 constexpr int num_eq = MSWEval::numWellEq;
2035 std::array<double, num_eq> tmp;
2036 for (
int ii = 0; ii < num_seg; ++ii) {
2037 const auto start = it + num_seg * num_eq;
2038 std::copy(start, start + num_eq, tmp.begin());
2039 this->primary_variables_.setValue(ii, tmp);
2041 return num_seg * num_eq;
Definition AquiferInterface.hpp:35
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:37
virtual void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition MultisegmentWell_impl.hpp:220
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
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:39