21#include <opm/common/Exceptions.hpp>
22#include <opm/common/OpmLog/OpmLog.hpp>
24#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
25#include <opm/input/eclipse/Units/Units.hpp>
27#include <opm/material/densead/EvaluationFormat.hpp>
29#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
30#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
31#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
36#if HAVE_CUDA || HAVE_OPENCL
37#include <opm/simulators/linalg/bda/WellContributions.hpp>
44 template <
typename TypeTag>
45 MultisegmentWell<TypeTag>::
46 MultisegmentWell(
const Well& well,
47 const ParallelWellInfo& pw_info,
49 const ModelParameters& param,
50 const RateConverterType& rate_converter,
51 const int pvtRegionIdx,
52 const int num_components,
54 const int index_of_well,
55 const std::vector<PerforationData>& perf_data)
56 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
57 , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
59 , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
62 if constexpr (has_solvent) {
63 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
66 if constexpr (has_polymer) {
67 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
70 if constexpr (Base::has_energy) {
71 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
74 if constexpr (Base::has_foam) {
75 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
78 if constexpr (Base::has_brine) {
79 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
82 if constexpr (Base::has_watVapor) {
83 OPM_THROW(std::runtime_error,
"water evaporation is not supported by multisegment well yet");
86 if(this->rsRvInj() > 0) {
87 OPM_THROW(std::runtime_error,
88 "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
89 " \n See (WCONINJE item 10 / WCONHIST item 8)");
91 if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
92 OPM_THROW(std::runtime_error,
"water + gas case not supported by multisegment well yet");
101 template <
typename TypeTag>
103 MultisegmentWell<TypeTag>::
104 init(
const PhaseUsage* phase_usage_arg,
105 const std::vector<double>& depth_arg,
106 const double gravity_arg,
108 const std::vector< Scalar >& B_avg,
109 const bool changed_to_open_this_step)
111 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
123 this->initMatrixAndVectors(num_cells);
126 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
127 const int cell_idx = this->well_cells_[perf];
128 this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
136 template <
typename TypeTag>
138 MultisegmentWell<TypeTag>::
139 initPrimaryVariablesEvaluation()
141 this->primary_variables_.init();
148 template <
typename TypeTag>
150 MultisegmentWell<TypeTag>::
151 updatePrimaryVariables(
const SummaryState& summary_state,
152 const WellState& well_state,
155 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
156 this->primary_variables_.update(well_state, stop_or_zero_rate_target);
164 template <
typename TypeTag>
175 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
176 this->segments_.perforations(),
178 this->scaleSegmentPressuresWithBhp(well_state);
185 template <
typename TypeTag>
189 const std::vector<double>&
B_avg,
193 return this->MSWEval::getWellConvergence(well_state,
196 this->param_.max_residual_allowed_,
197 this->param_.tolerance_wells_,
198 this->param_.relaxed_tolerance_flow_well_,
199 this->param_.tolerance_pressure_ms_wells_,
200 this->param_.relaxed_tolerance_pressure_ms_well_,
208 template <
typename TypeTag>
211 apply(
const BVector& x, BVector&
Ax)
const
213 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
217 if (this->param_.matrix_add_well_contributions_) {
229 template <
typename TypeTag>
234 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
243 template <
typename TypeTag>
251 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
256 this->linSys_.recoverSolutionWell(x,
xw);
264 template <
typename TypeTag>
269 std::vector<double>& well_potentials,
272 const int np = this->number_of_phases_;
273 well_potentials.resize(
np, 0.0);
276 if (this->wellIsStopped()) {
279 this->operability_status_.has_negative_potentials =
false;
284 const auto&
ws = well_state.well(this->index_of_well_);
285 if (this->isInjector()) {
286 const Well::InjectorCMode&
current =
ws.injection_cmode;
287 if (
current == Well::InjectorCMode::THP) {
290 if (
current == Well::InjectorCMode::BHP) {
294 const Well::ProducerCMode&
current =
ws.production_cmode;
295 if (
current == Well::ProducerCMode::THP) {
298 if (
current == Well::ProducerCMode::BHP) {
305 const double sign = this->isInjector() ? 1.0:-1.0;
306 for (
int phase = 0; phase <
np; ++phase){
313 for (
int phase = 0; phase <
np; ++phase){
314 well_potentials[phase] =
sign *
ws.surface_rates[phase];
320 debug_cost_counter_ = 0;
322 const auto& summaryState = ebosSimulator.vanguard().summaryState();
324 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials,
deferred_logger);
326 well_potentials = computeWellPotentialWithTHP(
329 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
330 + this->name() +
": " + std::to_string(debug_cost_counter_));
332 const double sign = this->isInjector() ? 1.0:-1.0;
334 for (
int phase = 0; phase <
np; ++phase){
335 well_potentials[phase] *=
sign;
340 this->operability_status_.has_negative_potentials =
true;
341 const std::string
msg = std::string(
"well ") + this->name() + std::string(
": has non negative potentials is not operable");
349 template<
typename TypeTag>
356 if (this->well_ecl_.isInjector()) {
357 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
360 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
361 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
365 template<
typename TypeTag>
367 MultisegmentWell<TypeTag>::
368 computeWellRatesWithBhp(
const Simulator& ebosSimulator,
370 std::vector<double>& well_flux,
371 DeferredLogger& deferred_logger)
const
374 const int np = this->number_of_phases_;
376 well_flux.resize(np, 0.0);
377 const bool allow_cf = this->getAllowCrossFlow();
378 const int nseg = this->numberOfSegments();
379 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
380 const auto& ws = well_state.well(this->indexOfWell());
381 auto segments_copy = ws.segments;
382 segments_copy.scale_pressure(bhp);
383 const auto& segment_pressure = segments_copy.pressure;
384 for (
int seg = 0; seg < nseg; ++seg) {
385 for (
const int perf : this->segments_.perforations()[seg]) {
386 const int cell_idx = this->well_cells_[perf];
387 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
389 std::vector<Scalar> mob(this->num_components_, 0.);
390 getMobilityScalar(ebosSimulator, perf, mob);
391 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
392 const double Tw = this->well_index_[perf] * trans_mult;
394 const Scalar seg_pressure = segment_pressure[seg];
395 std::vector<Scalar> cq_s(this->num_components_, 0.);
396 computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
397 allow_cf, cq_s, deferred_logger);
399 for(
int p = 0; p < np; ++p) {
400 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
404 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
408 template<
typename TypeTag>
410 MultisegmentWell<TypeTag>::
411 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
413 std::vector<double>& well_flux,
414 DeferredLogger& deferred_logger)
const
418 MultisegmentWell<TypeTag> well_copy(*
this);
419 well_copy.debug_cost_counter_ = 0;
422 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
423 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
424 auto& ws = well_state_copy.well(this->index_of_well_);
427 const auto& summary_state = ebosSimulator.vanguard().summaryState();
428 auto inj_controls = well_copy.well_ecl_.isInjector()
429 ? well_copy.well_ecl_.injectionControls(summary_state)
430 : Well::InjectionControls(0);
431 auto prod_controls = well_copy.well_ecl_.isProducer()
432 ? well_copy.well_ecl_.productionControls(summary_state) :
433 Well::ProductionControls(0);
436 if (well_copy.well_ecl_.isInjector()) {
437 inj_controls.bhp_limit = bhp;
438 ws.injection_cmode = Well::InjectorCMode::BHP;
440 prod_controls.bhp_limit = bhp;
441 ws.production_cmode = Well::ProducerCMode::BHP;
444 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
447 const int np = this->number_of_phases_;
449 for (
int phase = 0; phase < np; ++phase){
450 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
453 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
454 for (
int phase = 0; phase < np; ++phase) {
455 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
458 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
459 this->segments_.perforations(),
462 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
463 const double dt = ebosSimulator.timeStepSize();
465 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
470 well_flux.resize(np, 0.0);
471 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
472 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
473 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
475 debug_cost_counter_ += well_copy.debug_cost_counter_;
480 template<
typename TypeTag>
482 MultisegmentWell<TypeTag>::
483 computeWellPotentialWithTHP(
484 const WellState& well_state,
485 const Simulator& ebos_simulator,
486 DeferredLogger& deferred_logger)
const
488 std::vector<double> potentials(this->number_of_phases_, 0.0);
489 const auto& summary_state = ebos_simulator.vanguard().summaryState();
491 const auto& well = this->well_ecl_;
492 if (well.isInjector()){
493 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
494 if (bhp_at_thp_limit) {
495 const auto& controls = well.injectionControls(summary_state);
496 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
497 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
498 deferred_logger.debug(
"Converged thp based potential calculation for well "
499 + this->name() +
", at bhp = " + std::to_string(bhp));
501 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
502 "Failed in getting converged thp based potential calculation for well "
503 + this->name() +
". Instead the bhp based value is used");
504 const auto& controls = well.injectionControls(summary_state);
505 const double bhp = controls.bhp_limit;
506 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
509 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
510 well_state, ebos_simulator, summary_state, deferred_logger);
511 if (bhp_at_thp_limit) {
512 const auto& controls = well.productionControls(summary_state);
513 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
514 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
515 deferred_logger.debug(
"Converged thp based potential calculation for well "
516 + this->name() +
", at bhp = " + std::to_string(bhp));
518 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
519 "Failed in getting converged thp based potential calculation for well "
520 + this->name() +
". Instead the bhp based value is used");
521 const auto& controls = well.productionControls(summary_state);
522 const double bhp = controls.bhp_limit;
523 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
532 template <
typename TypeTag>
534 MultisegmentWell<TypeTag>::
535 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
537 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
539 std::vector<double> kr(this->number_of_phases_, 0.0);
540 std::vector<double> density(this->number_of_phases_, 0.0);
542 const int cell_idx = this->well_cells_[perf];
543 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
544 const auto& fs = intQuants.fluidState();
549 if (pu.phase_used[Water]) {
550 const int water_pos = pu.phase_pos[Water];
551 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
552 sum_kr += kr[water_pos];
553 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
556 if (pu.phase_used[Oil]) {
557 const int oil_pos = pu.phase_pos[Oil];
558 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
559 sum_kr += kr[oil_pos];
560 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
563 if (pu.phase_used[Gas]) {
564 const int gas_pos = pu.phase_pos[Gas];
565 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
566 sum_kr += kr[gas_pos];
567 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
570 assert(sum_kr != 0.);
573 double average_density = 0.;
574 for (
int p = 0; p < this->number_of_phases_; ++p) {
575 average_density += kr[p] * density[p];
577 average_density /= sum_kr;
579 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
587 template <
typename TypeTag>
589 MultisegmentWell<TypeTag>::
590 computeInitialSegmentFluids(
const Simulator& ebos_simulator)
592 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
594 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
595 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
596 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
605 template <
typename TypeTag>
607 MultisegmentWell<TypeTag>::
608 updateWellState(
const SummaryState& summary_state,
609 const BVectorWell& dwells,
610 WellState& well_state,
611 DeferredLogger& deferred_logger,
612 const double relaxation_factor)
614 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
616 const double dFLimit = this->param_.dwell_fraction_max_;
617 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
618 const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
619 this->primary_variables_.updateNewton(dwells,
622 stop_or_zero_rate_target,
623 max_pressure_change);
625 this->primary_variables_.copyToWellState(*
this, getRefDensity(), stop_or_zero_rate_target,
626 well_state, deferred_logger);
627 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
634 template <
typename TypeTag>
636 MultisegmentWell<TypeTag>::
637 calculateExplicitQuantities(
const Simulator& ebosSimulator,
638 const WellState& well_state,
639 DeferredLogger& deferred_logger)
641 const auto& summary_state = ebosSimulator.vanguard().summaryState();
642 updatePrimaryVariables(summary_state, well_state, deferred_logger);
643 initPrimaryVariablesEvaluation();
644 computePerfCellPressDiffs(ebosSimulator);
645 computeInitialSegmentFluids(ebosSimulator);
652 template<
typename TypeTag>
654 MultisegmentWell<TypeTag>::
655 updateProductivityIndex(
const Simulator& ebosSimulator,
656 const WellProdIndexCalculator& wellPICalc,
657 WellState& well_state,
658 DeferredLogger& deferred_logger)
const
660 auto fluidState = [&ebosSimulator,
this](
const int perf)
662 const auto cell_idx = this->well_cells_[perf];
663 return ebosSimulator.model()
664 .intensiveQuantities(cell_idx, 0).fluidState();
667 const int np = this->number_of_phases_;
668 auto setToZero = [np](
double* x) ->
void
670 std::fill_n(x, np, 0.0);
673 auto addVector = [np](
const double* src,
double* dest) ->
void
675 std::transform(src, src + np, dest, dest, std::plus<>{});
678 auto& ws = well_state.well(this->index_of_well_);
679 auto& perf_data = ws.perf_data;
680 auto* connPI = perf_data.prod_index.data();
681 auto* wellPI = ws.productivity_index.data();
685 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
686 auto subsetPerfID = 0;
688 for (
const auto& perf : *this->perf_data_){
689 auto allPerfID = perf.ecl_index;
691 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
693 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
696 std::vector<Scalar> mob(this->num_components_, 0.0);
697 getMobilityScalar(ebosSimulator,
static_cast<int>(subsetPerfID), mob);
699 const auto& fs = fluidState(subsetPerfID);
702 if (this->isInjector()) {
703 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
704 mob, connPI, deferred_logger);
707 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
710 addVector(connPI, wellPI);
716 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
717 "Internal logic error in processing connections for PI/II");
724 template<
typename TypeTag>
726 MultisegmentWell<TypeTag>::
727 addWellContributions(SparseMatrixAdapter& jacobian)
const
729 this->linSys_.extract(jacobian);
733 template<
typename TypeTag>
735 MultisegmentWell<TypeTag>::
736 addWellPressureEquations(PressureMatrix& jacobian,
737 const BVector& weights,
738 const int pressureVarIndex,
739 const bool use_well_weights,
740 const WellState& well_state)
const
743 this->linSys_.extractCPRPressureMatrix(jacobian,
753 template<
typename TypeTag>
754 template<
class Value>
756 MultisegmentWell<TypeTag>::
757 computePerfRate(
const Value& pressure_cell,
760 const std::vector<Value>& b_perfcells,
761 const std::vector<Value>& mob_perfcells,
764 const Value& segment_pressure,
765 const Value& segment_density,
766 const bool& allow_cf,
767 const std::vector<Value>& cmix_s,
768 std::vector<Value>& cq_s,
770 double& perf_dis_gas_rate,
771 double& perf_vap_oil_rate,
772 DeferredLogger& deferred_logger)
const
775 const Value perf_seg_press_diff = this->gravity() * segment_density *
776 this->segments_.perforation_depth_diff(perf);
778 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
782 perf_press = segment_pressure + perf_seg_press_diff;
785 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
788 const Value drawdown = cell_press_at_perf - perf_press;
791 if ( drawdown > 0.0) {
793 if (!allow_cf && this->isInjector()) {
798 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
799 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
800 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
803 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
804 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
805 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
806 const Value cq_s_oil = cq_s[oilCompIdx];
807 const Value cq_s_gas = cq_s[gasCompIdx];
808 cq_s[gasCompIdx] += rs * cq_s_oil;
809 cq_s[oilCompIdx] += rv * cq_s_gas;
813 if (!allow_cf && this->isProducer()) {
818 Value total_mob = mob_perfcells[0];
819 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
820 total_mob += mob_perfcells[comp_idx];
824 const Value cqt_i = - Tw * (total_mob * drawdown);
827 Value volume_ratio = 0.0;
828 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
829 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
830 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
833 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
834 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
835 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
840 const Value d = 1.0 - rv * rs;
842 if (getValue(d) == 0.0) {
843 OPM_DEFLOG_THROW(NumericalProblem,
844 fmt::format(
"Zero d value obtained for well {} "
845 "during flux calculation with rs {} and rv {}",
846 this->name(), rs, rv),
850 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
851 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
853 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
854 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
856 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
857 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
858 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
860 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
861 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
862 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
866 Value cqt_is = cqt_i / volume_ratio;
867 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
868 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
873 if (this->isProducer()) {
874 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
875 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
876 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
885 const double d = 1.0 - getValue(rv) * getValue(rs);
888 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
891 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
896 template <
typename TypeTag>
898 MultisegmentWell<TypeTag>::
899 computePerfRateEval(
const IntensiveQuantities& int_quants,
900 const std::vector<EvalWell>& mob_perfcells,
904 const EvalWell& segment_pressure,
905 const bool& allow_cf,
906 std::vector<EvalWell>& cq_s,
907 EvalWell& perf_press,
908 double& perf_dis_gas_rate,
909 double& perf_vap_oil_rate,
910 DeferredLogger& deferred_logger)
const
913 const auto& fs = int_quants.fluidState();
915 const EvalWell pressure_cell = this->extendEval(this->getPerfCellPressure(fs));
916 const EvalWell rs = this->extendEval(fs.Rs());
917 const EvalWell rv = this->extendEval(fs.Rv());
920 std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
922 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
923 if (!FluidSystem::phaseIsActive(phaseIdx)) {
927 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
928 b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
931 std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
932 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
933 cmix_s[comp_idx] = this->primary_variables_.surfaceVolumeFraction(seg, comp_idx);
936 this->computePerfRate(pressure_cell,
944 this->segments_.density(seg),
956 template <
typename TypeTag>
958 MultisegmentWell<TypeTag>::
959 computePerfRateScalar(
const IntensiveQuantities& int_quants,
960 const std::vector<Scalar>& mob_perfcells,
964 const Scalar& segment_pressure,
965 const bool& allow_cf,
966 std::vector<Scalar>& cq_s,
967 DeferredLogger& deferred_logger)
const
970 const auto& fs = int_quants.fluidState();
972 const Scalar pressure_cell = getValue(this->getPerfCellPressure(fs));
973 const Scalar rs = getValue(fs.Rs());
974 const Scalar rv = getValue(fs.Rv());
977 std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
979 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
980 if (!FluidSystem::phaseIsActive(phaseIdx)) {
984 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
985 b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
988 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
989 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
990 cmix_s[comp_idx] = getValue(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
993 Scalar perf_dis_gas_rate = 0.0;
994 Scalar perf_vap_oil_rate = 0.0;
995 Scalar perf_press = 0.0;
997 this->computePerfRate(pressure_cell,
1005 getValue(this->segments_.density(seg)),
1015 template <
typename TypeTag>
1017 MultisegmentWell<TypeTag>::
1018 computeSegmentFluidProperties(
const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
1027 EvalWell temperature;
1028 EvalWell saltConcentration;
1034 int pvt_region_index;
1037 const int cell_idx = this->well_cells_[0];
1038 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1039 const auto& fs = intQuants.fluidState();
1040 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1041 saltConcentration = this->extendEval(fs.saltConcentration());
1042 pvt_region_index = fs.pvtRegionIndex();
1045 this->segments_.computeFluidProperties(temperature,
1047 this->primary_variables_,
1056 template <
typename TypeTag>
1058 MultisegmentWell<TypeTag>::
1059 getMobilityEval(
const Simulator& ebosSimulator,
1061 std::vector<EvalWell>& mob)
const
1064 const int cell_idx = this->well_cells_[perf];
1065 assert (
int(mob.size()) == this->num_components_);
1066 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1067 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1071 const int satid = this->saturation_table_number_[perf] - 1;
1072 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1073 if( satid == satid_elem ) {
1075 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1076 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1080 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1081 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
1088 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1089 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
1090 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1093 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1096 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1097 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1101 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1102 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1108 template <
typename TypeTag>
1110 MultisegmentWell<TypeTag>::
1111 getMobilityScalar(
const Simulator& ebosSimulator,
1113 std::vector<Scalar>& mob)
const
1116 const int cell_idx = this->well_cells_[perf];
1117 assert (
int(mob.size()) == this->num_components_);
1118 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1119 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1123 const int satid = this->saturation_table_number_[perf] - 1;
1124 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1125 if( satid == satid_elem ) {
1127 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1128 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1132 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1133 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
1140 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1141 std::array<Scalar,3> relativePerms = { 0.0, 0.0, 0.0 };
1142 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1145 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1148 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1149 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1153 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1154 mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
1162 template<
typename TypeTag>
1164 MultisegmentWell<TypeTag>::
1165 getRefDensity()
const
1167 return this->segments_.getRefDensity();
1170 template<
typename TypeTag>
1172 MultisegmentWell<TypeTag>::
1173 checkOperabilityUnderBHPLimit(
const WellState& ,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1175 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1176 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1179 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1180 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1183 double total_ipr_mass_rate = 0.0;
1184 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1186 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1190 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1191 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1193 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1194 total_ipr_mass_rate += ipr_rate * rho;
1196 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1197 this->operability_status_.operable_under_only_bhp_limit =
false;
1201 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1205 std::vector<double> well_rates_bhp_limit;
1206 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1208 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1210 this->getRefDensity(),
1211 this->wellEcl().alq_value(),
1213 const double thp_limit = this->getTHPConstraint(summaryState);
1214 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1215 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1226 this->operability_status_.operable_under_only_bhp_limit =
true;
1227 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1233 template<
typename TypeTag>
1235 MultisegmentWell<TypeTag>::
1236 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1241 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1242 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1244 const int nseg = this->numberOfSegments();
1245 std::vector<double> seg_dp(nseg, 0.0);
1246 for (
int seg = 0; seg < nseg; ++seg) {
1248 const double dp = this->getSegmentDp(seg,
1249 this->segments_.density(seg).value(),
1252 for (
const int perf : this->segments_.perforations()[seg]) {
1253 std::vector<Scalar> mob(this->num_components_, 0.0);
1256 getMobilityScalar(ebos_simulator, perf, mob);
1258 const int cell_idx = this->well_cells_[perf];
1259 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1260 const auto& fs = int_quantities.fluidState();
1262 const double perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1264 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1265 const double pressure_cell = this->getPerfCellPressure(fs).value();
1268 std::vector<double> b_perf(this->num_components_);
1269 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1270 if (!FluidSystem::phaseIsActive(phase)) {
1273 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1274 b_perf[comp_idx] = fs.invB(phase).value();
1278 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1279 const double pressure_diff = pressure_cell - h_perf;
1282 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1283 deferred_logger.debug(
"CROSSFLOW_IPR",
1284 "cross flow found when updateIPR for well " + this->name());
1288 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1290 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1291 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1292 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1293 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1294 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1295 ipr_b_perf[comp_idx] += tw_mob;
1299 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1300 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1301 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1302 const double rs = (fs.Rs()).value();
1303 const double rv = (fs.Rv()).value();
1305 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1306 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1308 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1309 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1311 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1312 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1314 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1315 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1318 for (
size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1319 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1320 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1326 template<
typename TypeTag>
1328 MultisegmentWell<TypeTag>::
1329 checkOperabilityUnderTHPLimit(
1330 const Simulator& ebos_simulator,
1331 const WellState& well_state,
1332 DeferredLogger& deferred_logger)
1334 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1335 const auto obtain_bhp = this->isProducer()
1336 ? computeBhpAtThpLimitProd(
1337 well_state, ebos_simulator, summaryState, deferred_logger)
1338 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1341 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1343 const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1344 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1346 const double thp_limit = this->getTHPConstraint(summaryState);
1347 if (this->isProducer() && *obtain_bhp < thp_limit) {
1348 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1349 +
" bars is SMALLER than thp limit "
1350 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1351 +
" bars as a producer for well " + this->name();
1352 deferred_logger.debug(msg);
1354 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1355 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1356 +
" bars is LARGER than thp limit "
1357 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1358 +
" bars as a injector for well " + this->name();
1359 deferred_logger.debug(msg);
1364 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1365 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1366 if (!this->wellIsStopped()) {
1367 const double thp_limit = this->getTHPConstraint(summaryState);
1368 deferred_logger.debug(
" could not find bhp value at thp limit "
1369 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1370 +
" bar for well " + this->name() +
", the well might need to be closed ");
1379 template<
typename TypeTag>
1381 MultisegmentWell<TypeTag>::
1382 iterateWellEqWithControl(
const Simulator& ebosSimulator,
1384 const Well::InjectionControls& inj_controls,
1385 const Well::ProductionControls& prod_controls,
1386 WellState& well_state,
1387 const GroupState& group_state,
1388 DeferredLogger& deferred_logger)
1390 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1392 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1396 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1401 std::vector<std::vector<Scalar> > residual_history;
1402 std::vector<double> measure_history;
1405 double relaxation_factor = 1.;
1406 const double min_relaxation_factor = 0.6;
1407 bool converged =
false;
1408 int stagnate_count = 0;
1409 bool relax_convergence =
false;
1410 this->regularize_ =
false;
1411 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1413 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1415 const BVectorWell dx_well = this->linSys_.solve();
1417 if (it > this->param_.strict_inner_iter_wells_) {
1418 relax_convergence =
true;
1419 this->regularize_ =
true;
1422 const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
1423 if (report.converged()) {
1430 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1434 residual_history.push_back(residuals);
1435 measure_history.push_back(this->getResidualMeasureValue(well_state,
1436 residual_history[it],
1437 this->param_.tolerance_wells_,
1438 this->param_.tolerance_pressure_ms_wells_,
1443 bool is_oscillate =
false;
1444 bool is_stagnate =
false;
1446 this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
1450 if (is_oscillate || is_stagnate) {
1452 std::ostringstream sstr;
1453 if (relaxation_factor == min_relaxation_factor) {
1456 if (stagnate_count == 6) {
1457 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1458 const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger,
true);
1459 if (reportStag.converged()) {
1461 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1462 deferred_logger.debug(sstr.str());
1469 const double reduction_mutliplier = 0.9;
1470 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1474 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1478 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1480 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1482 this->regularize_ =
true;
1483 deferred_logger.debug(sstr.str());
1485 const auto& summary_state = ebosSimulator.vanguard().summaryState();
1486 updateWellState(summary_state, dx_well, well_state, deferred_logger, relaxation_factor);
1487 initPrimaryVariablesEvaluation();
1492 std::ostringstream sstr;
1493 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1494 if (relax_convergence)
1495 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1496 deferred_logger.debug(sstr.str());
1498 std::ostringstream sstr;
1499 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1500#define EXTRA_DEBUG_MSW 0
1502 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1503 for (
int i = 0; i < it; ++i) {
1504 const auto& residual = residual_history[i];
1505 sstr <<
" residual at " << i <<
"th iteration ";
1506 for (
const auto& res : residual) {
1509 sstr <<
" " << measure_history[i] <<
" \n";
1512#undef EXTRA_DEBUG_MSW
1513 deferred_logger.debug(sstr.str());
1523 template<
typename TypeTag>
1525 MultisegmentWell<TypeTag>::
1526 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
1528 const Well::InjectionControls& inj_controls,
1529 const Well::ProductionControls& prod_controls,
1530 WellState& well_state,
1531 const GroupState& group_state,
1532 DeferredLogger& deferred_logger)
1534 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1537 this->segments_.updateUpwindingSegments(this->primary_variables_);
1540 computeSegmentFluidProperties(ebosSimulator, deferred_logger);
1543 this->linSys_.clear();
1545 auto& ws = well_state.well(this->index_of_well_);
1546 ws.dissolved_gas_rate = 0;
1547 ws.dissolved_gas_rate_in_water = 0;
1548 ws.vaporized_oil_rate = 0;
1549 ws.vaporized_wat_rate = 0;
1556 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1558 const int nseg = this->numberOfSegments();
1560 for (
int seg = 0; seg < nseg; ++seg) {
1564 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1569 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1571 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1572 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1573 - segment_fluid_initial_[seg][comp_idx]) / dt;
1574 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1575 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1580 const int seg_upwind = this->segments_.upwinding_segment(seg);
1581 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1582 const EvalWell segment_rate =
1583 this->primary_variables_.getSegmentRateUpwinding(seg,
1586 this->well_efficiency_factor_;
1587 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1588 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1594 for (
const int inlet : this->segments_.inlets()[seg]) {
1595 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1596 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1597 const EvalWell inlet_rate =
1598 this->primary_variables_.getSegmentRateUpwinding(inlet,
1601 this->well_efficiency_factor_;
1602 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1603 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1609 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1610 auto& perf_data = ws.perf_data;
1611 auto& perf_rates = perf_data.phase_rates;
1612 auto& perf_press_state = perf_data.pressure;
1613 for (
const int perf : this->segments_.perforations()[seg]) {
1614 const int cell_idx = this->well_cells_[perf];
1615 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, 0);
1616 std::vector<EvalWell> mob(this->num_components_, 0.0);
1617 getMobilityEval(ebosSimulator, perf, mob);
1618 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1619 const double Tw = this->well_index_[perf] * trans_mult;
1620 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1621 EvalWell perf_press;
1622 double perf_dis_gas_rate = 0.;
1623 double perf_vap_oil_rate = 0.;
1624 computePerfRateEval(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
1627 if (this->isProducer()) {
1628 ws.dissolved_gas_rate += perf_dis_gas_rate;
1629 ws.vaporized_oil_rate += perf_vap_oil_rate;
1633 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1634 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1636 perf_press_state[perf] = perf_press.value();
1638 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1640 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1642 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1644 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1645 assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
1651 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1652 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1653 MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
1654 assembleControlEq(well_state,
1661 this->primary_variables_,
1665 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1666 this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
1670 this->linSys_.createSolver();
1676 template<
typename TypeTag>
1678 MultisegmentWell<TypeTag>::
1679 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1681 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1685 template<
typename TypeTag>
1687 MultisegmentWell<TypeTag>::
1688 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1690 bool all_drawdown_wrong_direction =
true;
1691 const int nseg = this->numberOfSegments();
1693 for (
int seg = 0; seg < nseg; ++seg) {
1694 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1695 for (
const int perf : this->segments_.perforations()[seg]) {
1697 const int cell_idx = this->well_cells_[perf];
1698 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1699 const auto& fs = intQuants.fluidState();
1702 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
1704 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1706 const double pressure_cell = this->getPerfCellPressure(fs).value();
1707 const double perf_press = pressure_cell - cell_perf_press_diff;
1710 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1715 if ( (drawdown < 0. && this->isInjector()) ||
1716 (drawdown > 0. && this->isProducer()) ) {
1717 all_drawdown_wrong_direction =
false;
1723 return all_drawdown_wrong_direction;
1729 template<
typename TypeTag>
1731 MultisegmentWell<TypeTag>::
1732 updateWaterThroughput(
const double , WellState& )
const
1740 template<
typename TypeTag>
1741 typename MultisegmentWell<TypeTag>::EvalWell
1742 MultisegmentWell<TypeTag>::
1743 getSegmentSurfaceVolume(
const Simulator& ebos_simulator,
const int seg_idx)
const
1745 EvalWell temperature;
1746 EvalWell saltConcentration;
1747 int pvt_region_index;
1751 const int cell_idx = this->well_cells_[0];
1752 const auto& intQuants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1753 const auto& fs = intQuants.fluidState();
1754 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1755 saltConcentration = this->extendEval(fs.saltConcentration());
1756 pvt_region_index = fs.pvtRegionIndex();
1759 return this->segments_.getSurfaceVolume(temperature,
1761 this->primary_variables_,
1767 template<
typename TypeTag>
1768 std::optional<double>
1769 MultisegmentWell<TypeTag>::
1770 computeBhpAtThpLimitProd(
const WellState& well_state,
1771 const Simulator& ebos_simulator,
1772 const SummaryState& summary_state,
1773 DeferredLogger& deferred_logger)
const
1775 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
1778 this->getALQ(well_state),
1784 template<
typename TypeTag>
1785 std::optional<double>
1786 MultisegmentWell<TypeTag>::
1787 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
1788 const SummaryState& summary_state,
1789 const double alq_value,
1790 DeferredLogger& deferred_logger)
const
1793 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1799 std::vector<double> rates(3);
1800 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1804 auto bhpAtLimit = WellBhpThpCalculator(*this).
1805 computeBhpAtThpLimitProd(frates,
1807 this->maxPerfPress(ebos_simulator),
1808 this->getRefDensity(),
1810 this->getTHPConstraint(summary_state),
1816 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1820 std::vector<double> rates(3);
1821 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1825 return WellBhpThpCalculator(*this).
1826 computeBhpAtThpLimitProd(fratesIter,
1828 this->maxPerfPress(ebos_simulator),
1829 this->getRefDensity(),
1831 this->getTHPConstraint(summary_state),
1835 template<
typename TypeTag>
1836 std::optional<double>
1837 MultisegmentWell<TypeTag>::
1838 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
1839 const SummaryState& summary_state,
1840 DeferredLogger& deferred_logger)
const
1843 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1849 std::vector<double> rates(3);
1850 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1854 auto bhpAtLimit = WellBhpThpCalculator(*this).
1855 computeBhpAtThpLimitInj(frates,
1857 this->getRefDensity(),
1866 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1870 std::vector<double> rates(3);
1871 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1875 return WellBhpThpCalculator(*this).
1876 computeBhpAtThpLimitInj(fratesIter,
1878 this->getRefDensity(),
1889 template<
typename TypeTag>
1891 MultisegmentWell<TypeTag>::
1892 maxPerfPress(
const Simulator& ebos_simulator)
const
1894 double max_pressure = 0.0;
1895 const int nseg = this->numberOfSegments();
1896 for (
int seg = 0; seg < nseg; ++seg) {
1897 for (
const int perf : this->segments_.perforations()[seg]) {
1898 const int cell_idx = this->well_cells_[perf];
1899 const auto& int_quants = ebos_simulator.model().intensiveQuantities(cell_idx, 0);
1900 const auto& fs = int_quants.fluidState();
1901 double pressure_cell = this->getPerfCellPressure(fs).value();
1902 max_pressure = std::max(max_pressure, pressure_cell);
1905 return max_pressure;
1912 template<
typename TypeTag>
1919 std::vector<Scalar>
well_q_s(this->num_components_, 0.0);
1920 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1921 const int nseg = this->numberOfSegments();
1925 for (
const int perf : this->segments_.perforations()[
seg]) {
1928 std::vector<Scalar>
mob(this->num_components_, 0.0);
1929 getMobilityScalar(ebosSimulator,
perf,
mob);
1932 std::vector<Scalar>
cq_s(this->num_components_, 0.0);
1934 for (
int comp = 0;
comp < this->num_components_; ++
comp) {
1946 template<
typename TypeTag>
1950 const std::function<
double(
const double)>&
connPICalc,
1951 const std::vector<Scalar>&
mobility,
1955 const int np = this->number_of_phases_;
1956 for (
int p = 0;
p <
np; ++
p) {
1960 mobility[ this->flowPhaseToEbosCompIdx(
p) ]
1961 *
fs.invB(this->flowPhaseToEbosPhaseIdx(
p)).value();
1966 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1967 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1969 const auto io = pu.phase_pos[Oil];
1970 const auto ig = pu.phase_pos[Gas];
1972 const auto vapoil = connPI[ig] * fs.Rv().value();
1973 const auto disgas = connPI[io] * fs.Rs().value();
1975 connPI[io] += vapoil;
1976 connPI[ig] += disgas;
1984 template<
typename TypeTag>
1986 MultisegmentWell<TypeTag>::
1987 computeConnLevelInjInd(
const typename MultisegmentWell<TypeTag>::FluidState& fs,
1988 const Phase preferred_phase,
1989 const std::function<
double(
const double)>& connIICalc,
1990 const std::vector<Scalar>& mobility,
1992 DeferredLogger& deferred_logger)
const
1998 if (preferred_phase == Phase::GAS) {
1999 phase_pos = pu.phase_pos[Gas];
2001 else if (preferred_phase == Phase::OIL) {
2002 phase_pos = pu.phase_pos[Oil];
2004 else if (preferred_phase == Phase::WATER) {
2005 phase_pos = pu.phase_pos[Water];
2008 OPM_DEFLOG_THROW(NotImplemented,
2009 fmt::format(
"Unsupported Injector Type ({}) "
2010 "for well {} during connection I.I. calculation",
2011 static_cast<int>(preferred_phase), this->name()),
2015 const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2016 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
2022 template<
typename TypeTag>
2024 MultisegmentWell<TypeTag>::
2025 updateWellStateWithTHPTargetProd(
const Simulator& ebos_simulator,
2026 WellState& well_state,
2027 DeferredLogger& deferred_logger)
const
2029 const auto& summary_state = ebos_simulator.vanguard().summaryState();
2031 auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
2032 ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
2033 if (bhp_at_thp_limit) {
2034 std::vector<double> rates(this->number_of_phases_, 0.0);
2035 computeWellRatesWithBhpIterations(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
2036 auto& ws = well_state.well(this->name());
2037 ws.surface_rates = rates;
2038 ws.bhp = *bhp_at_thp_limit;
2039 ws.thp = this->getTHPConstraint(summary_state);
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:211
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