22 #include <opm/simulators/wells/MSWellHelpers.hpp>
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
25 #include <opm/common/OpmLog/OpmLog.hpp>
30 #if HAVE_CUDA || HAVE_OPENCL
31 #include <opm/simulators/linalg/bda/WellContributions.hpp>
38 template <
typename TypeTag>
39 MultisegmentWell<TypeTag>::
40 MultisegmentWell(
const Well& well,
41 const ParallelWellInfo& pw_info,
43 const ModelParameters& param,
44 const RateConverterType& rate_converter,
45 const int pvtRegionIdx,
46 const int num_components,
48 const int index_of_well,
49 const std::vector<PerforationData>& perf_data)
50 : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
51 , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
52 , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
55 if constexpr (has_solvent) {
56 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
59 if constexpr (has_polymer) {
60 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
63 if constexpr (Base::has_energy) {
64 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
67 if constexpr (Base::has_foam) {
68 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
71 if constexpr (Base::has_brine) {
72 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
80 template <
typename TypeTag>
82 MultisegmentWell<TypeTag>::
83 init(
const PhaseUsage* phase_usage_arg,
84 const std::vector<double>& depth_arg,
85 const double gravity_arg,
87 const std::vector< Scalar >& B_avg)
89 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
101 this->initMatrixAndVectors(num_cells);
104 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
105 const int cell_idx = this->well_cells_[perf];
106 this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
114 template <
typename TypeTag>
116 MultisegmentWell<TypeTag>::
117 initPrimaryVariablesEvaluation()
const
119 this->MSWEval::initPrimaryVariablesEvaluation();
126 template <
typename TypeTag>
128 MultisegmentWell<TypeTag>::
129 updatePrimaryVariables(
const WellState& well_state, DeferredLogger& )
const
131 this->MSWEval::updatePrimaryVariables(well_state);
139 template <
typename TypeTag>
147 Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
150 this->scaleSegmentRatesWithWellRates(well_state);
151 this->scaleSegmentPressuresWithBhp(well_state);
158 template <
typename TypeTag>
162 const std::vector<double>& B_avg,
164 const bool relax_tolerance)
const
166 return this->MSWEval::getWellConvergence(well_state,
169 this->param_.max_residual_allowed_,
170 this->param_.tolerance_wells_,
171 this->param_.relaxed_tolerance_flow_well_,
172 this->param_.tolerance_pressure_ms_wells_,
173 this->param_.relaxed_tolerance_pressure_ms_well_,
181 template <
typename TypeTag>
184 apply(
const BVector& x, BVector& Ax)
const
186 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
188 if ( this->param_.matrix_add_well_contributions_ )
193 BVectorWell Bx(this->duneB_.N());
195 this->duneB_.mv(x, Bx);
198 const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx);
201 this->duneC_.mmtv(invDBx,Ax);
208 template <
typename TypeTag>
211 apply(BVector& r)
const
213 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
216 const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
218 this->duneC_.mmtv(invDrw, r);
223 template <
typename TypeTag>
230 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
233 this->recoverSolutionWell(x, xw);
234 updateWellState(xw, well_state, deferred_logger);
241 template <
typename TypeTag>
246 std::vector<double>& well_potentials,
249 const int np = this->number_of_phases_;
250 well_potentials.resize(np, 0.0);
253 if (this->wellIsStopped()) {
258 bool thp_controlled_well =
false;
259 bool bhp_controlled_well =
false;
260 const auto& ws = well_state.well(this->index_of_well_);
261 if (this->isInjector()) {
262 const Well::InjectorCMode& current = ws.injection_cmode;
263 if (current == Well::InjectorCMode::THP) {
264 thp_controlled_well =
true;
266 if (current == Well::InjectorCMode::BHP) {
267 bhp_controlled_well =
true;
270 const Well::ProducerCMode& current = ws.production_cmode;
271 if (current == Well::ProducerCMode::THP) {
272 thp_controlled_well =
true;
274 if (current == Well::ProducerCMode::BHP) {
275 bhp_controlled_well =
true;
278 if (thp_controlled_well || bhp_controlled_well) {
280 double total_rate = 0.0;
281 for (
int phase = 0; phase < np; ++phase){
282 total_rate += ws.surface_rates[phase];
287 if (std::abs(total_rate) > 0) {
288 for (
int phase = 0; phase < np; ++phase){
289 well_potentials[phase] = ws.surface_rates[phase];
295 debug_cost_counter_ = 0;
297 const auto& summaryState = ebosSimulator.vanguard().summaryState();
298 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
299 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
301 well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
303 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
304 + this->name() +
": " + std::to_string(debug_cost_counter_));
310 template<
typename TypeTag>
314 std::vector<double>& well_flux,
317 if (this->well_ecl_.isInjector()) {
318 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
319 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
321 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
322 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
326 template<
typename TypeTag>
328 MultisegmentWell<TypeTag>::
329 computeWellRatesWithBhp(
const Simulator& ebosSimulator,
331 std::vector<double>& well_flux,
332 DeferredLogger& deferred_logger)
const
335 const int np = this->number_of_phases_;
337 well_flux.resize(np, 0.0);
338 const bool allow_cf = this->getAllowCrossFlow();
339 const int nseg = this->numberOfSegments();
340 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
341 const auto& ws = well_state.well(this->indexOfWell());
342 auto segments_copy = ws.segments;
343 segments_copy.scale_pressure(bhp);
344 const auto& segment_pressure = segments_copy.pressure;
345 for (
int seg = 0; seg < nseg; ++seg) {
346 for (
const int perf : this->segment_perforations_[seg]) {
347 const int cell_idx = this->well_cells_[perf];
348 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
350 std::vector<Scalar> mob(this->num_components_, 0.);
351 getMobilityScalar(ebosSimulator, perf, mob);
352 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
353 const double Tw = this->well_index_[perf] * trans_mult;
355 const Scalar seg_pressure = segment_pressure[seg];
356 std::vector<Scalar> cq_s(this->num_components_, 0.);
357 computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
358 allow_cf, cq_s, 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_;
409 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
410 for (
int phase = 0; phase < np; ++phase){
411 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
413 well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
415 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
416 const double dt = ebosSimulator.timeStepSize();
418 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
423 well_flux.resize(np, 0.0);
424 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
425 const EvalWell rate = well_copy.getQs(compIdx);
426 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
428 debug_cost_counter_ += well_copy.debug_cost_counter_;
433 template<
typename TypeTag>
435 MultisegmentWell<TypeTag>::
436 computeWellPotentialWithTHP(
const Simulator& ebos_simulator,
437 DeferredLogger& deferred_logger)
const
439 std::vector<double> potentials(this->number_of_phases_, 0.0);
440 const auto& summary_state = ebos_simulator.vanguard().summaryState();
442 const auto& well = this->well_ecl_;
443 if (well.isInjector()){
444 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
445 if (bhp_at_thp_limit) {
446 const auto& controls = well.injectionControls(summary_state);
447 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
448 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
449 deferred_logger.debug(
"Converged thp based potential calculation for well "
450 + this->name() +
", at bhp = " + std::to_string(bhp));
452 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
453 "Failed in getting converged thp based potential calculation for well "
454 + this->name() +
". Instead the bhp based value is used");
455 const auto& controls = well.injectionControls(summary_state);
456 const double bhp = controls.bhp_limit;
457 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
460 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
461 if (bhp_at_thp_limit) {
462 const auto& controls = well.productionControls(summary_state);
463 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
464 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
465 deferred_logger.debug(
"Converged thp based potential calculation for well "
466 + this->name() +
", at bhp = " + std::to_string(bhp));
468 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
469 "Failed in getting converged thp based potential calculation for well "
470 + this->name() +
". Instead the bhp based value is used");
471 const auto& controls = well.productionControls(summary_state);
472 const double bhp = controls.bhp_limit;
473 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
482 template <
typename TypeTag>
484 MultisegmentWell<TypeTag>::
485 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
487 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
491 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
493 updateWellState(dx_well, well_state, deferred_logger);
500 template <
typename TypeTag>
502 MultisegmentWell<TypeTag>::
503 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
505 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
507 std::vector<double> kr(this->number_of_phases_, 0.0);
508 std::vector<double> density(this->number_of_phases_, 0.0);
510 const int cell_idx = this->well_cells_[perf];
511 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
512 const auto& fs = intQuants.fluidState();
517 if (pu.phase_used[Water]) {
518 const int water_pos = pu.phase_pos[Water];
519 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
520 sum_kr += kr[water_pos];
521 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
524 if (pu.phase_used[Oil]) {
525 const int oil_pos = pu.phase_pos[Oil];
526 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
527 sum_kr += kr[oil_pos];
528 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
531 if (pu.phase_used[Gas]) {
532 const int gas_pos = pu.phase_pos[Gas];
533 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
534 sum_kr += kr[gas_pos];
535 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
538 assert(sum_kr != 0.);
541 double average_density = 0.;
542 for (
int p = 0; p < this->number_of_phases_; ++p) {
543 average_density += kr[p] * density[p];
545 average_density /= sum_kr;
547 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
555 template <
typename TypeTag>
557 MultisegmentWell<TypeTag>::
558 computeInitialSegmentFluids(
const Simulator& ebos_simulator)
560 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
562 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
563 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
564 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
573 template <
typename TypeTag>
575 MultisegmentWell<TypeTag>::
576 updateWellState(
const BVectorWell& dwells,
577 WellState& well_state,
578 DeferredLogger& deferred_logger,
579 const double relaxation_factor)
const
581 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
583 const double dFLimit = this->param_.dwell_fraction_max_;
584 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
585 this->MSWEval::updateWellState(dwells,
588 max_pressure_change);
590 this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
591 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
598 template <
typename TypeTag>
600 MultisegmentWell<TypeTag>::
601 calculateExplicitQuantities(
const Simulator& ebosSimulator,
602 const WellState& well_state,
603 DeferredLogger& deferred_logger)
605 updatePrimaryVariables(well_state, deferred_logger);
606 initPrimaryVariablesEvaluation();
607 computePerfCellPressDiffs(ebosSimulator);
608 computeInitialSegmentFluids(ebosSimulator);
615 template<
typename TypeTag>
617 MultisegmentWell<TypeTag>::
618 updateProductivityIndex(
const Simulator& ebosSimulator,
619 const WellProdIndexCalculator& wellPICalc,
620 WellState& well_state,
621 DeferredLogger& deferred_logger)
const
623 auto fluidState = [&ebosSimulator,
this](
const int perf)
625 const auto cell_idx = this->well_cells_[perf];
626 return ebosSimulator.model()
627 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
630 const int np = this->number_of_phases_;
631 auto setToZero = [np](
double* x) ->
void
633 std::fill_n(x, np, 0.0);
636 auto addVector = [np](
const double* src,
double* dest) ->
void
638 std::transform(src, src + np, dest, dest, std::plus<>{});
641 auto& ws = well_state.well(this->index_of_well_);
642 auto& perf_data = ws.perf_data;
643 auto* connPI = perf_data.prod_index.data();
644 auto* wellPI = ws.productivity_index.data();
648 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
649 auto subsetPerfID = 0;
651 for (
const auto& perf : *this->perf_data_){
652 auto allPerfID = perf.ecl_index;
654 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
656 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
659 std::vector<Scalar> mob(this->num_components_, 0.0);
660 getMobilityScalar(ebosSimulator,
static_cast<int>(subsetPerfID), mob);
662 const auto& fs = fluidState(subsetPerfID);
665 if (this->isInjector()) {
666 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
667 mob, connPI, deferred_logger);
670 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
673 addVector(connPI, wellPI);
679 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
680 "Internal logic error in processing connections for PI/II");
687 template<
typename TypeTag>
689 MultisegmentWell<TypeTag>::
690 addWellContributions(SparseMatrixAdapter& jacobian)
const
692 const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(this->duneD_, this->duneDSolver_);
702 for (
size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
703 for (
auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
704 const auto row_index = colC.index();
705 for (
size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
706 for (
auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
707 const auto col_index = colB.index();
708 OffDiagMatrixBlockWellType tmp1;
709 Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
710 typename SparseMatrixAdapter::MatrixBlock tmp2;
711 Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
712 jacobian.addToBlock(row_index, col_index, tmp2);
721 template<
typename TypeTag>
722 template<
class Value>
724 MultisegmentWell<TypeTag>::
725 computePerfRate(
const Value& pressure_cell,
728 const std::vector<Value>& b_perfcells,
729 const std::vector<Value>& mob_perfcells,
732 const Value& segment_pressure,
733 const Value& segment_density,
734 const bool& allow_cf,
735 const std::vector<Value>& cmix_s,
736 std::vector<Value>& cq_s,
738 double& perf_dis_gas_rate,
739 double& perf_vap_oil_rate,
740 DeferredLogger& deferred_logger)
const
743 const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
745 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
747 perf_press = pressure_cell - cell_perf_press_diff;
750 const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
753 if ( drawdown > 0.0) {
755 if (!allow_cf && this->isInjector()) {
760 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
761 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
762 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
765 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
766 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
767 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
768 const Value cq_s_oil = cq_s[oilCompIdx];
769 const Value cq_s_gas = cq_s[gasCompIdx];
770 cq_s[gasCompIdx] += rs * cq_s_oil;
771 cq_s[oilCompIdx] += rv * cq_s_gas;
775 if (!allow_cf && this->isProducer()) {
780 Value total_mob = mob_perfcells[0];
781 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
782 total_mob += mob_perfcells[comp_idx];
786 const Value cqt_i = - Tw * (total_mob * drawdown);
789 Value volume_ratio = 0.0;
790 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
791 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
792 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
795 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
796 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
797 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
802 const Value d = 1.0 - rv * rs;
804 if (getValue(d) == 0.0) {
805 OPM_DEFLOG_THROW(NumericalIssue,
"Zero d value obtained for well " << this->name()
806 <<
" during flux calculation"
807 <<
" with rs " << rs <<
" and rv " << rv, deferred_logger);
810 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
811 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
813 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
814 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
816 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
817 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
818 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
820 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
821 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
822 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
826 Value cqt_is = cqt_i / volume_ratio;
827 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
828 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
833 if (this->isProducer()) {
834 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
835 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
836 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
845 const double d = 1.0 - getValue(rv) * getValue(rs);
848 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
851 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
856 template <
typename TypeTag>
858 MultisegmentWell<TypeTag>::
859 computePerfRateEval(
const IntensiveQuantities& int_quants,
860 const std::vector<EvalWell>& mob_perfcells,
864 const EvalWell& segment_pressure,
865 const bool& allow_cf,
866 std::vector<EvalWell>& cq_s,
867 EvalWell& perf_press,
868 double& perf_dis_gas_rate,
869 double& perf_vap_oil_rate,
870 DeferredLogger& deferred_logger)
const
873 const auto& fs = int_quants.fluidState();
875 const EvalWell pressure_cell = this->extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
876 const EvalWell rs = this->extendEval(fs.Rs());
877 const EvalWell rv = this->extendEval(fs.Rv());
880 std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
882 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
883 if (!FluidSystem::phaseIsActive(phaseIdx)) {
887 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
888 b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
891 std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
892 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
893 cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
896 this->computePerfRate(pressure_cell,
904 this->segment_densities_[seg],
916 template <
typename TypeTag>
918 MultisegmentWell<TypeTag>::
919 computePerfRateScalar(
const IntensiveQuantities& int_quants,
920 const std::vector<Scalar>& mob_perfcells,
924 const Scalar& segment_pressure,
925 const bool& allow_cf,
926 std::vector<Scalar>& cq_s,
927 DeferredLogger& deferred_logger)
const
930 const auto& fs = int_quants.fluidState();
932 const Scalar pressure_cell = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
933 const Scalar rs = getValue(fs.Rs());
934 const Scalar rv = getValue(fs.Rv());
937 std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
939 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
940 if (!FluidSystem::phaseIsActive(phaseIdx)) {
944 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
945 b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
948 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
949 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
950 cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
953 Scalar perf_dis_gas_rate = 0.0;
954 Scalar perf_vap_oil_rate = 0.0;
955 Scalar perf_press = 0.0;
957 this->computePerfRate(pressure_cell,
965 getValue(this->segment_densities_[seg]),
975 template <
typename TypeTag>
977 MultisegmentWell<TypeTag>::
978 computeSegmentFluidProperties(
const Simulator& ebosSimulator)
987 EvalWell temperature;
988 EvalWell saltConcentration;
994 int pvt_region_index;
997 const int cell_idx = this->well_cells_[0];
998 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
999 const auto& fs = intQuants.fluidState();
1000 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1001 saltConcentration = this->extendEval(fs.saltConcentration());
1002 pvt_region_index = fs.pvtRegionIndex();
1005 this->MSWEval::computeSegmentFluidProperties(temperature,
1014 template <
typename TypeTag>
1016 MultisegmentWell<TypeTag>::
1017 getMobilityEval(
const Simulator& ebosSimulator,
1019 std::vector<EvalWell>& mob)
const
1022 const int cell_idx = this->well_cells_[perf];
1023 assert (
int(mob.size()) == this->num_components_);
1024 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1025 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1029 const int satid = this->saturation_table_number_[perf] - 1;
1030 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1031 if( satid == satid_elem ) {
1033 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1034 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1038 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1039 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
1046 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1047 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
1048 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1051 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1054 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1055 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1059 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1060 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1066 template <
typename TypeTag>
1068 MultisegmentWell<TypeTag>::
1069 getMobilityScalar(
const Simulator& ebosSimulator,
1071 std::vector<Scalar>& mob)
const
1074 const int cell_idx = this->well_cells_[perf];
1075 assert (
int(mob.size()) == this->num_components_);
1076 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1077 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1081 const int satid = this->saturation_table_number_[perf] - 1;
1082 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1083 if( satid == satid_elem ) {
1085 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1086 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1090 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1091 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
1098 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1099 Scalar relativePerms[3] = { 0.0, 0.0, 0.0 };
1100 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1103 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1106 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1107 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1111 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1112 mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
1120 template<
typename TypeTag>
1122 MultisegmentWell<TypeTag>::
1123 getRefDensity()
const
1125 return this->segment_densities_[0].value();
1128 template<
typename TypeTag>
1130 MultisegmentWell<TypeTag>::
1131 checkOperabilityUnderBHPLimitProducer(
const WellState& ,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1133 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1134 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1137 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1138 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1143 for (
int p = 0; p < this->number_of_phases_; ++p) {
1144 temp += this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
1147 this->operability_status_.operable_under_only_bhp_limit =
false;
1151 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1155 std::vector<double> well_rates_bhp_limit;
1156 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1158 const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
1160 const double thp_limit = this->getTHPConstraint(summaryState);
1162 if (thp < thp_limit) {
1163 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1174 this->operability_status_.operable_under_only_bhp_limit =
true;
1175 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1181 template<
typename TypeTag>
1183 MultisegmentWell<TypeTag>::
1184 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1190 if (this->isInjector()) {
1195 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1196 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1198 const int nseg = this->numberOfSegments();
1199 double seg_bhp_press_diff = 0;
1200 double ref_depth = this->ref_depth_;
1201 for (
int seg = 0; seg < nseg; ++seg) {
1203 const double segment_depth = this->segmentSet()[seg].depth();
1204 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
1205 ref_depth = segment_depth;
1206 seg_bhp_press_diff += dp;
1207 for (
const int perf : this->segment_perforations_[seg]) {
1208 std::vector<Scalar> mob(this->num_components_, 0.0);
1211 getMobilityScalar(ebos_simulator, perf, mob);
1213 const int cell_idx = this->well_cells_[perf];
1214 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1215 const auto& fs = int_quantities.fluidState();
1218 const double perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
1220 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1221 const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
1224 std::vector<double> b_perf(this->num_components_);
1225 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1226 if (!FluidSystem::phaseIsActive(phase)) {
1229 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1230 b_perf[comp_idx] = fs.invB(phase).value();
1234 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + seg_bhp_press_diff;
1235 const double pressure_diff = pressure_cell - h_perf;
1240 if (pressure_diff <= 0.) {
1241 deferred_logger.warning(
"NON_POSITIVE_DRAWDOWN_IPR",
1242 "non-positive drawdown found when updateIPR for well " + this->name());
1246 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1251 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1252 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1253 for (
int p = 0; p < this->number_of_phases_; ++p) {
1254 const double tw_mob = tw_perf * mob[p] * b_perf[p];
1255 ipr_a_perf[p] += tw_mob * pressure_diff;
1256 ipr_b_perf[p] += tw_mob;
1260 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1261 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1262 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1263 const double rs = (fs.Rs()).value();
1264 const double rv = (fs.Rv()).value();
1266 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1267 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1269 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1270 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1272 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1273 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1275 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1276 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1279 for (
int p = 0; p < this->number_of_phases_; ++p) {
1281 this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
1282 this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
1288 template<
typename TypeTag>
1290 MultisegmentWell<TypeTag>::
1291 checkOperabilityUnderTHPLimitProducer(
const Simulator& ebos_simulator,
const WellState& , DeferredLogger& deferred_logger)
1293 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1294 const auto obtain_bhp = computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
1297 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1299 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1300 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1302 const double thp_limit = this->getTHPConstraint(summaryState);
1303 if (*obtain_bhp < thp_limit) {
1304 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1305 +
" bars is SMALLER than thp limit "
1306 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1307 +
" bars as a producer for well " + this->name();
1308 deferred_logger.debug(msg);
1313 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1314 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1315 if (!this->wellIsStopped()) {
1316 const double thp_limit = this->getTHPConstraint(summaryState);
1317 deferred_logger.debug(
" could not find bhp value at thp limit "
1318 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1319 +
" bar for well " + this->name() +
", the well might need to be closed ");
1328 template<
typename TypeTag>
1330 MultisegmentWell<TypeTag>::
1331 iterateWellEqWithControl(
const Simulator& ebosSimulator,
1333 const Well::InjectionControls& inj_controls,
1334 const Well::ProductionControls& prod_controls,
1335 WellState& well_state,
1336 const GroupState& group_state,
1337 DeferredLogger& deferred_logger)
1339 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1341 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1342 const WellState well_state0 = well_state;
1343 const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
1344 std::vector<std::vector<Scalar> > residual_history;
1345 std::vector<double> measure_history;
1348 double relaxation_factor = 1.;
1349 const double min_relaxation_factor = 0.6;
1350 bool converged =
false;
1351 int stagnate_count = 0;
1352 bool relax_convergence =
false;
1353 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1355 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1357 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
1359 if (it > this->param_.strict_inner_iter_wells_)
1360 relax_convergence =
true;
1362 const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
1363 if (report.converged()) {
1368 residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
1369 measure_history.push_back(this->getResidualMeasureValue(well_state,
1370 residual_history[it],
1371 this->param_.tolerance_wells_,
1372 this->param_.tolerance_pressure_ms_wells_,
1375 bool is_oscillate =
false;
1376 bool is_stagnate =
false;
1378 this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
1382 if (is_oscillate || is_stagnate) {
1384 std::ostringstream sstr;
1385 if (relaxation_factor == min_relaxation_factor) {
1388 if (stagnate_count == 6) {
1389 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1390 const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger,
true);
1391 if (reportStag.converged()) {
1393 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1394 deferred_logger.debug(sstr.str());
1401 const double reduction_mutliplier = 0.9;
1402 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1406 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1410 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1412 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1413 deferred_logger.debug(sstr.str());
1415 updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
1416 initPrimaryVariablesEvaluation();
1421 std::ostringstream sstr;
1422 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1423 if (relax_convergence)
1424 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1425 deferred_logger.debug(sstr.str());
1427 std::ostringstream sstr;
1428 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1429 #define EXTRA_DEBUG_MSW 0
1431 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1432 for (
int i = 0; i < it; ++i) {
1433 const auto& residual = residual_history[i];
1434 sstr <<
" residual at " << i <<
"th iteration ";
1435 for (
const auto& res : residual) {
1438 sstr <<
" " << measure_history[i] <<
" \n";
1441 deferred_logger.debug(sstr.str());
1451 template<
typename TypeTag>
1453 MultisegmentWell<TypeTag>::
1454 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
1456 const Well::InjectionControls& inj_controls,
1457 const Well::ProductionControls& prod_controls,
1458 WellState& well_state,
1459 const GroupState& group_state,
1460 DeferredLogger& deferred_logger)
1463 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1466 this->updateUpwindingSegments();
1469 computeSegmentFluidProperties(ebosSimulator);
1476 this->resWell_ = 0.0;
1478 this->duneDSolver_.reset();
1480 auto& ws = well_state.well(this->index_of_well_);
1481 ws.dissolved_gas_rate = 0;
1482 ws.vaporized_oil_rate = 0;
1489 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1491 const int nseg = this->numberOfSegments();
1493 for (
int seg = 0; seg < nseg; ++seg) {
1497 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1502 const Scalar regularization_factor = this->param_.regularization_factor_ms_wells_;
1504 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1505 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
1506 - segment_fluid_initial_[seg][comp_idx]) / dt;
1508 this->resWell_[seg][comp_idx] += accumulation_term.value();
1509 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1510 this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
1516 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1517 const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
1519 const int seg_upwind = this->upwinding_segments_[seg];
1522 this->resWell_[seg][comp_idx] -= segment_rate.value();
1523 this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq);
1524 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1525 this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
1527 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1528 this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
1536 for (
const int inlet : this->segment_inlets_[seg]) {
1537 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1538 const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
1540 const int inlet_upwind = this->upwinding_segments_[inlet];
1543 this->resWell_[seg][comp_idx] += inlet_rate.value();
1544 this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq);
1545 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1546 this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
1548 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1549 this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
1557 const EvalWell seg_pressure = this->getSegmentPressure(seg);
1558 auto& perf_data = ws.perf_data;
1559 auto& perf_rates = perf_data.phase_rates;
1560 auto& perf_press_state = perf_data.pressure;
1561 for (
const int perf : this->segment_perforations_[seg]) {
1562 const int cell_idx = this->well_cells_[perf];
1563 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1564 std::vector<EvalWell> mob(this->num_components_, 0.0);
1565 getMobilityEval(ebosSimulator, perf, mob);
1566 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1567 const double Tw = this->well_index_[perf] * trans_mult;
1568 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1569 EvalWell perf_press;
1570 double perf_dis_gas_rate = 0.;
1571 double perf_vap_oil_rate = 0.;
1572 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);
1575 if (this->isProducer()) {
1576 ws.dissolved_gas_rate += perf_dis_gas_rate;
1577 ws.vaporized_oil_rate += perf_vap_oil_rate;
1581 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1582 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1584 perf_press_state[perf] = perf_press.value();
1586 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1588 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1590 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1593 this->resWell_[seg][comp_idx] += cq_s_effective.value();
1596 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1599 this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq);
1602 this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
1605 for (
int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
1607 this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
1614 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1615 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1616 this->assembleControlEq(well_state,
1625 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1626 this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
1634 template<
typename TypeTag>
1636 MultisegmentWell<TypeTag>::
1637 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1639 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1643 template<
typename TypeTag>
1645 MultisegmentWell<TypeTag>::
1646 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1648 bool all_drawdown_wrong_direction =
true;
1649 const int nseg = this->numberOfSegments();
1651 for (
int seg = 0; seg < nseg; ++seg) {
1652 const EvalWell segment_pressure = this->getSegmentPressure(seg);
1653 for (
const int perf : this->segment_perforations_[seg]) {
1655 const int cell_idx = this->well_cells_[perf];
1656 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1657 const auto& fs = intQuants.fluidState();
1660 const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
1662 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1664 const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
1665 const double perf_press = pressure_cell - cell_perf_press_diff;
1668 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1673 if ( (drawdown < 0. && this->isInjector()) ||
1674 (drawdown > 0. && this->isProducer()) ) {
1675 all_drawdown_wrong_direction =
false;
1681 return all_drawdown_wrong_direction;
1687 template<
typename TypeTag>
1689 MultisegmentWell<TypeTag>::
1690 updateWaterThroughput(
const double , WellState& )
const
1698 template<
typename TypeTag>
1699 typename MultisegmentWell<TypeTag>::EvalWell
1700 MultisegmentWell<TypeTag>::
1701 getSegmentSurfaceVolume(
const Simulator& ebos_simulator,
const int seg_idx)
const
1703 EvalWell temperature;
1704 EvalWell saltConcentration;
1705 int pvt_region_index;
1709 const int cell_idx = this->well_cells_[0];
1710 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1711 const auto& fs = intQuants.fluidState();
1712 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1713 saltConcentration = this->extendEval(fs.saltConcentration());
1714 pvt_region_index = fs.pvtRegionIndex();
1717 return this->MSWEval::getSegmentSurfaceVolume(temperature,
1727 template<
typename TypeTag>
1728 std::optional<double>
1729 MultisegmentWell<TypeTag>::
1730 computeBhpAtThpLimitProd(
const Simulator& ebos_simulator,
1731 const SummaryState& summary_state,
1732 DeferredLogger& deferred_logger)
const
1735 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1741 std::vector<double> rates(3);
1742 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1746 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1747 computeBhpAtThpLimitProd(frates,
1749 maxPerfPress(ebos_simulator),
1756 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1760 std::vector<double> rates(3);
1761 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1765 return this->MultisegmentWellGeneric<Scalar>::
1766 computeBhpAtThpLimitProd(fratesIter,
1768 maxPerfPress(ebos_simulator),
1777 template<
typename TypeTag>
1778 std::optional<double>
1779 MultisegmentWell<TypeTag>::
1780 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
1781 const SummaryState& summary_state,
1782 DeferredLogger& deferred_logger)
const
1785 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1791 std::vector<double> rates(3);
1792 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1796 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1797 computeBhpAtThpLimitInj(frates,
1805 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1809 std::vector<double> rates(3);
1810 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1814 return this->MultisegmentWellGeneric<Scalar>::
1815 computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
1822 template<
typename TypeTag>
1824 MultisegmentWell<TypeTag>::
1825 maxPerfPress(
const Simulator& ebos_simulator)
const
1827 double max_pressure = 0.0;
1828 const int nseg = this->numberOfSegments();
1829 for (
int seg = 0; seg < nseg; ++seg) {
1830 for (
const int perf : this->segment_perforations_[seg]) {
1831 const int cell_idx = this->well_cells_[perf];
1832 const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1833 const auto& fs = int_quants.fluidState();
1834 double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
1835 max_pressure = std::max(max_pressure, pressure_cell);
1838 return max_pressure;
1845 template<
typename TypeTag>
1852 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
1853 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1854 const int nseg = this->numberOfSegments();
1855 for (
int seg = 0; seg < nseg; ++seg) {
1857 const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
1858 for (
const int perf : this->segment_perforations_[seg]) {
1859 const int cell_idx = this->well_cells_[perf];
1860 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1861 std::vector<Scalar> mob(this->num_components_, 0.0);
1862 getMobilityScalar(ebosSimulator, perf, mob);
1863 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1864 const double Tw = this->well_index_[perf] * trans_mult;
1865 std::vector<Scalar> cq_s(this->num_components_, 0.0);
1866 computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger);
1867 for (
int comp = 0; comp < this->num_components_; ++comp) {
1868 well_q_s[comp] += cq_s[comp];
1879 template<
typename TypeTag>
1883 const std::function<
double(
const double)>& connPICalc,
1884 const std::vector<Scalar>& mobility,
1885 double* connPI)
const
1888 const int np = this->number_of_phases_;
1889 for (
int p = 0; p < np; ++p) {
1892 const auto connMob =
1893 mobility[ this->flowPhaseToEbosCompIdx(p) ]
1894 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
1896 connPI[p] = connPICalc(connMob);
1899 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1900 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1902 const auto io = pu.phase_pos[Oil];
1903 const auto ig = pu.phase_pos[Gas];
1905 const auto vapoil = connPI[ig] * fs.Rv().value();
1906 const auto disgas = connPI[io] * fs.Rs().value();
1908 connPI[io] += vapoil;
1909 connPI[ig] += disgas;
1917 template<
typename TypeTag>
1919 MultisegmentWell<TypeTag>::
1920 computeConnLevelInjInd(
const typename MultisegmentWell<TypeTag>::FluidState& fs,
1921 const Phase preferred_phase,
1922 const std::function<
double(
const double)>& connIICalc,
1923 const std::vector<Scalar>& mobility,
1925 DeferredLogger& deferred_logger)
const
1931 if (preferred_phase == Phase::GAS) {
1932 phase_pos = pu.phase_pos[Gas];
1934 else if (preferred_phase == Phase::OIL) {
1935 phase_pos = pu.phase_pos[Oil];
1937 else if (preferred_phase == Phase::WATER) {
1938 phase_pos = pu.phase_pos[Water];
1941 OPM_DEFLOG_THROW(NotImplemented,
1942 "Unsupported Injector Type ("
1943 <<
static_cast<int>(preferred_phase)
1944 <<
") for well " << this->name()
1945 <<
" during connection I.I. calculation",
1949 const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1950 connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition: ConvergenceReport.hpp:36
Definition: DeferredLogger.hpp:57
Definition: GroupState.hpp:34
Definition: MultisegmentWell.hpp:39
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33