22 #include <opm/simulators/wells/MSWellHelpers.hpp>
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/input/eclipse/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()) {
256 this->operability_status_.has_negative_potentials =
false;
259 bool thp_controlled_well =
false;
260 bool bhp_controlled_well =
false;
261 const auto& ws = well_state.well(this->index_of_well_);
262 if (this->isInjector()) {
263 const Well::InjectorCMode& current = ws.injection_cmode;
264 if (current == Well::InjectorCMode::THP) {
265 thp_controlled_well =
true;
267 if (current == Well::InjectorCMode::BHP) {
268 bhp_controlled_well =
true;
271 const Well::ProducerCMode& current = ws.production_cmode;
272 if (current == Well::ProducerCMode::THP) {
273 thp_controlled_well =
true;
275 if (current == Well::ProducerCMode::BHP) {
276 bhp_controlled_well =
true;
279 if (thp_controlled_well || bhp_controlled_well) {
281 double total_rate = 0.0;
282 const double sign = this->isInjector() ? 1.0:-1.0;
283 for (
int phase = 0; phase < np; ++phase){
284 total_rate += sign * ws.surface_rates[phase];
289 if (total_rate > 0) {
290 for (
int phase = 0; phase < np; ++phase){
291 well_potentials[phase] = sign * ws.surface_rates[phase];
297 debug_cost_counter_ = 0;
299 const auto& summaryState = ebosSimulator.vanguard().summaryState();
300 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
301 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
303 well_potentials = computeWellPotentialWithTHP(
304 well_state, ebosSimulator, deferred_logger);
306 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
307 + this->name() +
": " + std::to_string(debug_cost_counter_));
309 const double sign = this->isInjector() ? 1.0:-1.0;
310 double total_potential = 0.0;
311 for (
int phase = 0; phase < np; ++phase){
312 well_potentials[phase] *= sign;
313 total_potential += well_potentials[phase];
315 if (total_potential < 0.0 && this->param_.check_well_operability_) {
317 this->operability_status_.has_negative_potentials =
true;
318 const std::string msg = std::string(
"well ") + this->name() + std::string(
": has non negative potentials is not operable");
319 deferred_logger.warning(
"NEGATIVE_POTENTIALS_INOPERABLE", msg);
326 template<
typename TypeTag>
330 std::vector<double>& well_flux,
333 if (this->well_ecl_.isInjector()) {
334 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
335 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
337 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
338 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
342 template<
typename TypeTag>
344 MultisegmentWell<TypeTag>::
345 computeWellRatesWithBhp(
const Simulator& ebosSimulator,
347 std::vector<double>& well_flux,
348 DeferredLogger& deferred_logger)
const
351 const int np = this->number_of_phases_;
353 well_flux.resize(np, 0.0);
354 const bool allow_cf = this->getAllowCrossFlow();
355 const int nseg = this->numberOfSegments();
356 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
357 const auto& ws = well_state.well(this->indexOfWell());
358 auto segments_copy = ws.segments;
359 segments_copy.scale_pressure(bhp);
360 const auto& segment_pressure = segments_copy.pressure;
361 for (
int seg = 0; seg < nseg; ++seg) {
362 for (
const int perf : this->segment_perforations_[seg]) {
363 const int cell_idx = this->well_cells_[perf];
364 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
366 std::vector<Scalar> mob(this->num_components_, 0.);
367 getMobilityScalar(ebosSimulator, perf, mob);
368 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
369 const double Tw = this->well_index_[perf] * trans_mult;
371 const Scalar seg_pressure = segment_pressure[seg];
372 std::vector<Scalar> cq_s(this->num_components_, 0.);
373 computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
374 allow_cf, cq_s, deferred_logger);
376 for(
int p = 0; p < np; ++p) {
377 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
381 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
385 template<
typename TypeTag>
387 MultisegmentWell<TypeTag>::
388 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
390 std::vector<double>& well_flux,
391 DeferredLogger& deferred_logger)
const
395 MultisegmentWell<TypeTag> well_copy(*
this);
396 well_copy.debug_cost_counter_ = 0;
399 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
400 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
401 auto& ws = well_state_copy.well(this->index_of_well_);
404 const auto& summary_state = ebosSimulator.vanguard().summaryState();
405 auto inj_controls = well_copy.well_ecl_.isInjector()
406 ? well_copy.well_ecl_.injectionControls(summary_state)
407 : Well::InjectionControls(0);
408 auto prod_controls = well_copy.well_ecl_.isProducer()
409 ? well_copy.well_ecl_.productionControls(summary_state) :
410 Well::ProductionControls(0);
413 if (well_copy.well_ecl_.isInjector()) {
414 inj_controls.bhp_limit = bhp;
415 ws.injection_cmode = Well::InjectorCMode::BHP;
417 prod_controls.bhp_limit = bhp;
418 ws.production_cmode = Well::ProducerCMode::BHP;
421 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
424 const int np = this->number_of_phases_;
426 for (
int phase = 0; phase < np; ++phase){
427 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
430 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
431 for (
int phase = 0; phase < np; ++phase) {
432 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
435 well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
437 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
438 const double dt = ebosSimulator.timeStepSize();
440 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
445 well_flux.resize(np, 0.0);
446 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
447 const EvalWell rate = well_copy.getQs(compIdx);
448 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
450 debug_cost_counter_ += well_copy.debug_cost_counter_;
455 template<
typename TypeTag>
457 MultisegmentWell<TypeTag>::
458 computeWellPotentialWithTHP(
459 const WellState& well_state,
460 const Simulator& ebos_simulator,
461 DeferredLogger& deferred_logger)
const
463 std::vector<double> potentials(this->number_of_phases_, 0.0);
464 const auto& summary_state = ebos_simulator.vanguard().summaryState();
466 const auto& well = this->well_ecl_;
467 if (well.isInjector()){
468 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
469 if (bhp_at_thp_limit) {
470 const auto& controls = well.injectionControls(summary_state);
471 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
472 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
473 deferred_logger.debug(
"Converged thp based potential calculation for well "
474 + this->name() +
", at bhp = " + std::to_string(bhp));
476 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
477 "Failed in getting converged thp based potential calculation for well "
478 + this->name() +
". Instead the bhp based value is used");
479 const auto& controls = well.injectionControls(summary_state);
480 const double bhp = controls.bhp_limit;
481 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
484 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
485 well_state, ebos_simulator, summary_state, deferred_logger);
486 if (bhp_at_thp_limit) {
487 const auto& controls = well.productionControls(summary_state);
488 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
489 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
490 deferred_logger.debug(
"Converged thp based potential calculation for well "
491 + this->name() +
", at bhp = " + std::to_string(bhp));
493 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
494 "Failed in getting converged thp based potential calculation for well "
495 + this->name() +
". Instead the bhp based value is used");
496 const auto& controls = well.productionControls(summary_state);
497 const double bhp = controls.bhp_limit;
498 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
507 template <
typename TypeTag>
509 MultisegmentWell<TypeTag>::
510 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
512 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
516 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
518 updateWellState(dx_well, well_state, deferred_logger);
525 template <
typename TypeTag>
527 MultisegmentWell<TypeTag>::
528 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
530 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
532 std::vector<double> kr(this->number_of_phases_, 0.0);
533 std::vector<double> density(this->number_of_phases_, 0.0);
535 const int cell_idx = this->well_cells_[perf];
536 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
537 const auto& fs = intQuants.fluidState();
542 if (pu.phase_used[Water]) {
543 const int water_pos = pu.phase_pos[Water];
544 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
545 sum_kr += kr[water_pos];
546 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
549 if (pu.phase_used[Oil]) {
550 const int oil_pos = pu.phase_pos[Oil];
551 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
552 sum_kr += kr[oil_pos];
553 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
556 if (pu.phase_used[Gas]) {
557 const int gas_pos = pu.phase_pos[Gas];
558 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
559 sum_kr += kr[gas_pos];
560 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
563 assert(sum_kr != 0.);
566 double average_density = 0.;
567 for (
int p = 0; p < this->number_of_phases_; ++p) {
568 average_density += kr[p] * density[p];
570 average_density /= sum_kr;
572 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
580 template <
typename TypeTag>
582 MultisegmentWell<TypeTag>::
583 computeInitialSegmentFluids(
const Simulator& ebos_simulator)
585 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
587 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
588 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
589 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
598 template <
typename TypeTag>
600 MultisegmentWell<TypeTag>::
601 updateWellState(
const BVectorWell& dwells,
602 WellState& well_state,
603 DeferredLogger& deferred_logger,
604 const double relaxation_factor)
const
606 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
608 const double dFLimit = this->param_.dwell_fraction_max_;
609 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
610 this->MSWEval::updateWellState(dwells,
613 max_pressure_change);
615 this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
616 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
623 template <
typename TypeTag>
625 MultisegmentWell<TypeTag>::
626 calculateExplicitQuantities(
const Simulator& ebosSimulator,
627 const WellState& well_state,
628 DeferredLogger& deferred_logger)
630 updatePrimaryVariables(well_state, deferred_logger);
631 initPrimaryVariablesEvaluation();
632 computePerfCellPressDiffs(ebosSimulator);
633 computeInitialSegmentFluids(ebosSimulator);
640 template<
typename TypeTag>
642 MultisegmentWell<TypeTag>::
643 updateProductivityIndex(
const Simulator& ebosSimulator,
644 const WellProdIndexCalculator& wellPICalc,
645 WellState& well_state,
646 DeferredLogger& deferred_logger)
const
648 auto fluidState = [&ebosSimulator,
this](
const int perf)
650 const auto cell_idx = this->well_cells_[perf];
651 return ebosSimulator.model()
652 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
655 const int np = this->number_of_phases_;
656 auto setToZero = [np](
double* x) ->
void
658 std::fill_n(x, np, 0.0);
661 auto addVector = [np](
const double* src,
double* dest) ->
void
663 std::transform(src, src + np, dest, dest, std::plus<>{});
666 auto& ws = well_state.well(this->index_of_well_);
667 auto& perf_data = ws.perf_data;
668 auto* connPI = perf_data.prod_index.data();
669 auto* wellPI = ws.productivity_index.data();
673 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
674 auto subsetPerfID = 0;
676 for (
const auto& perf : *this->perf_data_){
677 auto allPerfID = perf.ecl_index;
679 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
681 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
684 std::vector<Scalar> mob(this->num_components_, 0.0);
685 getMobilityScalar(ebosSimulator,
static_cast<int>(subsetPerfID), mob);
687 const auto& fs = fluidState(subsetPerfID);
690 if (this->isInjector()) {
691 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
692 mob, connPI, deferred_logger);
695 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
698 addVector(connPI, wellPI);
704 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
705 "Internal logic error in processing connections for PI/II");
712 template<
typename TypeTag>
714 MultisegmentWell<TypeTag>::
715 addWellContributions(SparseMatrixAdapter& jacobian)
const
717 const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(this->duneD_, this->duneDSolver_);
727 for (
size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
728 for (
auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
729 const auto row_index = colC.index();
730 for (
size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
731 for (
auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
732 const auto col_index = colB.index();
733 OffDiagMatrixBlockWellType tmp1;
734 Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
735 typename SparseMatrixAdapter::MatrixBlock tmp2;
736 Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
737 jacobian.addToBlock(row_index, col_index, tmp2);
746 template<
typename TypeTag>
747 template<
class Value>
749 MultisegmentWell<TypeTag>::
750 computePerfRate(
const Value& pressure_cell,
753 const std::vector<Value>& b_perfcells,
754 const std::vector<Value>& mob_perfcells,
757 const Value& segment_pressure,
758 const Value& segment_density,
759 const bool& allow_cf,
760 const std::vector<Value>& cmix_s,
761 std::vector<Value>& cq_s,
763 double& perf_dis_gas_rate,
764 double& perf_vap_oil_rate,
765 DeferredLogger& deferred_logger)
const
768 const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
770 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
772 perf_press = pressure_cell - cell_perf_press_diff;
775 const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
778 if ( drawdown > 0.0) {
780 if (!allow_cf && this->isInjector()) {
785 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
786 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
787 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
790 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
791 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
792 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
793 const Value cq_s_oil = cq_s[oilCompIdx];
794 const Value cq_s_gas = cq_s[gasCompIdx];
795 cq_s[gasCompIdx] += rs * cq_s_oil;
796 cq_s[oilCompIdx] += rv * cq_s_gas;
800 if (!allow_cf && this->isProducer()) {
805 Value total_mob = mob_perfcells[0];
806 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
807 total_mob += mob_perfcells[comp_idx];
811 const Value cqt_i = - Tw * (total_mob * drawdown);
814 Value volume_ratio = 0.0;
815 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
816 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
817 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
820 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
821 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
822 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
827 const Value d = 1.0 - rv * rs;
829 if (getValue(d) == 0.0) {
830 OPM_DEFLOG_THROW(NumericalIssue,
"Zero d value obtained for well " << this->name()
831 <<
" during flux calculation"
832 <<
" with rs " << rs <<
" and rv " << rv, deferred_logger);
835 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
836 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
838 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
839 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
841 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
842 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
843 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
845 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
846 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
847 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
851 Value cqt_is = cqt_i / volume_ratio;
852 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
853 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
858 if (this->isProducer()) {
859 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
860 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
861 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
870 const double d = 1.0 - getValue(rv) * getValue(rs);
873 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
876 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
881 template <
typename TypeTag>
883 MultisegmentWell<TypeTag>::
884 computePerfRateEval(
const IntensiveQuantities& int_quants,
885 const std::vector<EvalWell>& mob_perfcells,
889 const EvalWell& segment_pressure,
890 const bool& allow_cf,
891 std::vector<EvalWell>& cq_s,
892 EvalWell& perf_press,
893 double& perf_dis_gas_rate,
894 double& perf_vap_oil_rate,
895 DeferredLogger& deferred_logger)
const
898 const auto& fs = int_quants.fluidState();
900 const EvalWell pressure_cell = this->extendEval(this->getPerfCellPressure(fs));
901 const EvalWell rs = this->extendEval(fs.Rs());
902 const EvalWell rv = this->extendEval(fs.Rv());
905 std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
907 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
908 if (!FluidSystem::phaseIsActive(phaseIdx)) {
912 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
913 b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
916 std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
917 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
918 cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
921 this->computePerfRate(pressure_cell,
929 this->segment_densities_[seg],
941 template <
typename TypeTag>
943 MultisegmentWell<TypeTag>::
944 computePerfRateScalar(
const IntensiveQuantities& int_quants,
945 const std::vector<Scalar>& mob_perfcells,
949 const Scalar& segment_pressure,
950 const bool& allow_cf,
951 std::vector<Scalar>& cq_s,
952 DeferredLogger& deferred_logger)
const
955 const auto& fs = int_quants.fluidState();
957 const Scalar pressure_cell = getValue(this->getPerfCellPressure(fs));
958 const Scalar rs = getValue(fs.Rs());
959 const Scalar rv = getValue(fs.Rv());
962 std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
964 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
965 if (!FluidSystem::phaseIsActive(phaseIdx)) {
969 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
970 b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
973 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
974 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
975 cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
978 Scalar perf_dis_gas_rate = 0.0;
979 Scalar perf_vap_oil_rate = 0.0;
980 Scalar perf_press = 0.0;
982 this->computePerfRate(pressure_cell,
990 getValue(this->segment_densities_[seg]),
1000 template <
typename TypeTag>
1002 MultisegmentWell<TypeTag>::
1003 computeSegmentFluidProperties(
const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
1012 EvalWell temperature;
1013 EvalWell saltConcentration;
1019 int pvt_region_index;
1022 const int cell_idx = this->well_cells_[0];
1023 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1024 const auto& fs = intQuants.fluidState();
1025 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1026 saltConcentration = this->extendEval(fs.saltConcentration());
1027 pvt_region_index = fs.pvtRegionIndex();
1030 this->MSWEval::computeSegmentFluidProperties(temperature,
1040 template <
typename TypeTag>
1042 MultisegmentWell<TypeTag>::
1043 getMobilityEval(
const Simulator& ebosSimulator,
1045 std::vector<EvalWell>& mob)
const
1048 const int cell_idx = this->well_cells_[perf];
1049 assert (
int(mob.size()) == this->num_components_);
1050 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1051 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1055 const int satid = this->saturation_table_number_[perf] - 1;
1056 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1057 if( satid == satid_elem ) {
1059 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1060 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1064 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1065 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
1072 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1073 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
1074 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1077 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1080 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1081 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1085 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1086 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1092 template <
typename TypeTag>
1094 MultisegmentWell<TypeTag>::
1095 getMobilityScalar(
const Simulator& ebosSimulator,
1097 std::vector<Scalar>& mob)
const
1100 const int cell_idx = this->well_cells_[perf];
1101 assert (
int(mob.size()) == this->num_components_);
1102 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1103 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1107 const int satid = this->saturation_table_number_[perf] - 1;
1108 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1109 if( satid == satid_elem ) {
1111 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1112 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1116 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1117 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
1124 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1125 Scalar relativePerms[3] = { 0.0, 0.0, 0.0 };
1126 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1129 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1132 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1133 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1137 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1138 mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
1146 template<
typename TypeTag>
1148 MultisegmentWell<TypeTag>::
1149 getRefDensity()
const
1151 return this->segment_densities_[0].value();
1154 template<
typename TypeTag>
1156 MultisegmentWell<TypeTag>::
1157 checkOperabilityUnderBHPLimit(
const WellState& ,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1159 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1160 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1163 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1164 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1168 double ipr_rate = 0;
1169 for (
int p = 0; p < this->number_of_phases_; ++p) {
1170 ipr_rate += this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
1172 if ( (this->isProducer() && ipr_rate < 0.) || (this->isInjector() && ipr_rate > 0.) ) {
1173 this->operability_status_.operable_under_only_bhp_limit =
false;
1177 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1181 std::vector<double> well_rates_bhp_limit;
1182 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1184 const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
1185 const double thp_limit = this->getTHPConstraint(summaryState);
1186 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1187 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1198 this->operability_status_.operable_under_only_bhp_limit =
true;
1199 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1205 template<
typename TypeTag>
1207 MultisegmentWell<TypeTag>::
1208 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1213 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1214 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1216 const int nseg = this->numberOfSegments();
1217 double seg_bhp_press_diff = 0;
1218 double ref_depth = this->ref_depth_;
1219 for (
int seg = 0; seg < nseg; ++seg) {
1221 const double segment_depth = this->segmentSet()[seg].depth();
1222 const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
1223 ref_depth = segment_depth;
1224 seg_bhp_press_diff += dp;
1225 for (
const int perf : this->segment_perforations_[seg]) {
1226 std::vector<Scalar> mob(this->num_components_, 0.0);
1229 getMobilityScalar(ebos_simulator, perf, mob);
1231 const int cell_idx = this->well_cells_[perf];
1232 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1233 const auto& fs = int_quantities.fluidState();
1236 const double perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
1238 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1239 const double pressure_cell = this->getPerfCellPressure(fs).value();
1242 std::vector<double> b_perf(this->num_components_);
1243 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1244 if (!FluidSystem::phaseIsActive(phase)) {
1247 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1248 b_perf[comp_idx] = fs.invB(phase).value();
1252 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + seg_bhp_press_diff;
1253 const double pressure_diff = pressure_cell - h_perf;
1256 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1257 deferred_logger.debug(
"CROSSFLOW_IPR",
1258 "cross flow found when updateIPR for well " + this->name());
1262 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1267 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1268 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1269 for (
int p = 0; p < this->number_of_phases_; ++p) {
1270 const double tw_mob = tw_perf * mob[p] * b_perf[p];
1271 ipr_a_perf[p] += tw_mob * pressure_diff;
1272 ipr_b_perf[p] += tw_mob;
1276 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1277 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1278 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1279 const double rs = (fs.Rs()).value();
1280 const double rv = (fs.Rv()).value();
1282 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1283 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1285 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1286 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1288 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1289 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1291 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1292 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1295 for (
int p = 0; p < this->number_of_phases_; ++p) {
1297 this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
1298 this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
1304 template<
typename TypeTag>
1306 MultisegmentWell<TypeTag>::
1307 checkOperabilityUnderTHPLimit(
1308 const Simulator& ebos_simulator,
1309 const WellState& well_state,
1310 DeferredLogger& deferred_logger)
1312 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1313 const auto obtain_bhp = this->isProducer()
1314 ? computeBhpAtThpLimitProd(
1315 well_state, ebos_simulator, summaryState, deferred_logger)
1316 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1319 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1321 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1322 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1324 const double thp_limit = this->getTHPConstraint(summaryState);
1325 if (this->isProducer() && *obtain_bhp < thp_limit) {
1326 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1327 +
" bars is SMALLER than thp limit "
1328 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1329 +
" bars as a producer for well " + this->name();
1330 deferred_logger.debug(msg);
1332 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1333 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1334 +
" bars is LARGER than thp limit "
1335 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1336 +
" bars as a injector for well " + this->name();
1337 deferred_logger.debug(msg);
1342 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1343 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1344 if (!this->wellIsStopped()) {
1345 const double thp_limit = this->getTHPConstraint(summaryState);
1346 deferred_logger.debug(
" could not find bhp value at thp limit "
1347 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1348 +
" bar for well " + this->name() +
", the well might need to be closed ");
1357 template<
typename TypeTag>
1359 MultisegmentWell<TypeTag>::
1360 iterateWellEqWithControl(
const Simulator& ebosSimulator,
1362 const Well::InjectionControls& inj_controls,
1363 const Well::ProductionControls& prod_controls,
1364 WellState& well_state,
1365 const GroupState& group_state,
1366 DeferredLogger& deferred_logger)
1368 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1370 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1371 const WellState well_state0 = well_state;
1375 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1380 std::vector<std::vector<Scalar> > residual_history;
1381 std::vector<double> measure_history;
1384 double relaxation_factor = 1.;
1385 const double min_relaxation_factor = 0.6;
1386 bool converged =
false;
1387 int stagnate_count = 0;
1388 bool relax_convergence =
false;
1389 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1391 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1393 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
1395 if (it > this->param_.strict_inner_iter_wells_)
1396 relax_convergence =
true;
1398 const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
1399 if (report.converged()) {
1406 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1410 residual_history.push_back(residuals);
1411 measure_history.push_back(this->getResidualMeasureValue(well_state,
1412 residual_history[it],
1413 this->param_.tolerance_wells_,
1414 this->param_.tolerance_pressure_ms_wells_,
1419 bool is_oscillate =
false;
1420 bool is_stagnate =
false;
1422 this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
1426 if (is_oscillate || is_stagnate) {
1428 std::ostringstream sstr;
1429 if (relaxation_factor == min_relaxation_factor) {
1432 if (stagnate_count == 6) {
1433 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1434 const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger,
true);
1435 if (reportStag.converged()) {
1437 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1438 deferred_logger.debug(sstr.str());
1445 const double reduction_mutliplier = 0.9;
1446 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1450 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1454 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1456 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1457 deferred_logger.debug(sstr.str());
1459 updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
1460 initPrimaryVariablesEvaluation();
1465 std::ostringstream sstr;
1466 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1467 if (relax_convergence)
1468 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1469 deferred_logger.debug(sstr.str());
1471 std::ostringstream sstr;
1472 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1473 #define EXTRA_DEBUG_MSW 0
1475 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1476 for (
int i = 0; i < it; ++i) {
1477 const auto& residual = residual_history[i];
1478 sstr <<
" residual at " << i <<
"th iteration ";
1479 for (
const auto& res : residual) {
1482 sstr <<
" " << measure_history[i] <<
" \n";
1485 deferred_logger.debug(sstr.str());
1495 template<
typename TypeTag>
1497 MultisegmentWell<TypeTag>::
1498 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
1500 const Well::InjectionControls& inj_controls,
1501 const Well::ProductionControls& prod_controls,
1502 WellState& well_state,
1503 const GroupState& group_state,
1504 DeferredLogger& deferred_logger)
1507 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1510 this->updateUpwindingSegments();
1513 computeSegmentFluidProperties(ebosSimulator, deferred_logger);
1520 this->resWell_ = 0.0;
1522 this->duneDSolver_.reset();
1524 auto& ws = well_state.well(this->index_of_well_);
1525 ws.dissolved_gas_rate = 0;
1526 ws.vaporized_oil_rate = 0;
1533 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1535 const int nseg = this->numberOfSegments();
1537 for (
int seg = 0; seg < nseg; ++seg) {
1541 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1546 const Scalar regularization_factor = this->param_.regularization_factor_ms_wells_;
1548 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1549 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
1550 - segment_fluid_initial_[seg][comp_idx]) / dt;
1552 this->resWell_[seg][comp_idx] += accumulation_term.value();
1553 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1554 this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
1560 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1561 const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
1563 const int seg_upwind = this->upwinding_segments_[seg];
1566 this->resWell_[seg][comp_idx] -= segment_rate.value();
1567 this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq);
1568 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1569 this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
1571 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1572 this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
1580 for (
const int inlet : this->segment_inlets_[seg]) {
1581 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1582 const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
1584 const int inlet_upwind = this->upwinding_segments_[inlet];
1587 this->resWell_[seg][comp_idx] += inlet_rate.value();
1588 this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq);
1589 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1590 this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
1592 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1593 this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
1601 const EvalWell seg_pressure = this->getSegmentPressure(seg);
1602 auto& perf_data = ws.perf_data;
1603 auto& perf_rates = perf_data.phase_rates;
1604 auto& perf_press_state = perf_data.pressure;
1605 for (
const int perf : this->segment_perforations_[seg]) {
1606 const int cell_idx = this->well_cells_[perf];
1607 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1608 std::vector<EvalWell> mob(this->num_components_, 0.0);
1609 getMobilityEval(ebosSimulator, perf, mob);
1610 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1611 const double Tw = this->well_index_[perf] * trans_mult;
1612 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1613 EvalWell perf_press;
1614 double perf_dis_gas_rate = 0.;
1615 double perf_vap_oil_rate = 0.;
1616 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);
1619 if (this->isProducer()) {
1620 ws.dissolved_gas_rate += perf_dis_gas_rate;
1621 ws.vaporized_oil_rate += perf_vap_oil_rate;
1625 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1626 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1628 perf_press_state[perf] = perf_press.value();
1630 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1632 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1634 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1637 this->resWell_[seg][comp_idx] += cq_s_effective.value();
1640 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1643 this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq);
1646 this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
1649 for (
int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
1651 this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
1658 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1659 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1660 this->assembleControlEq(well_state,
1669 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1670 this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
1678 template<
typename TypeTag>
1680 MultisegmentWell<TypeTag>::
1681 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1683 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1687 template<
typename TypeTag>
1689 MultisegmentWell<TypeTag>::
1690 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1692 bool all_drawdown_wrong_direction =
true;
1693 const int nseg = this->numberOfSegments();
1695 for (
int seg = 0; seg < nseg; ++seg) {
1696 const EvalWell segment_pressure = this->getSegmentPressure(seg);
1697 for (
const int perf : this->segment_perforations_[seg]) {
1699 const int cell_idx = this->well_cells_[perf];
1700 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1701 const auto& fs = intQuants.fluidState();
1704 const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
1706 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1708 const double pressure_cell = this->getPerfCellPressure(fs).value();
1709 const double perf_press = pressure_cell - cell_perf_press_diff;
1712 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1717 if ( (drawdown < 0. && this->isInjector()) ||
1718 (drawdown > 0. && this->isProducer()) ) {
1719 all_drawdown_wrong_direction =
false;
1725 return all_drawdown_wrong_direction;
1731 template<
typename TypeTag>
1733 MultisegmentWell<TypeTag>::
1734 updateWaterThroughput(
const double , WellState& )
const
1742 template<
typename TypeTag>
1743 typename MultisegmentWell<TypeTag>::EvalWell
1744 MultisegmentWell<TypeTag>::
1745 getSegmentSurfaceVolume(
const Simulator& ebos_simulator,
const int seg_idx)
const
1747 EvalWell temperature;
1748 EvalWell saltConcentration;
1749 int pvt_region_index;
1753 const int cell_idx = this->well_cells_[0];
1754 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1755 const auto& fs = intQuants.fluidState();
1756 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1757 saltConcentration = this->extendEval(fs.saltConcentration());
1758 pvt_region_index = fs.pvtRegionIndex();
1761 return this->MSWEval::getSegmentSurfaceVolume(temperature,
1768 template<
typename TypeTag>
1769 std::optional<double>
1770 MultisegmentWell<TypeTag>::
1771 computeBhpAtThpLimitProd(
const WellState& well_state,
1772 const Simulator& ebos_simulator,
1773 const SummaryState& summary_state,
1774 DeferredLogger& deferred_logger)
const
1776 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
1780 this->getALQ(well_state));
1785 template<
typename TypeTag>
1786 std::optional<double>
1787 MultisegmentWell<TypeTag>::
1788 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
1789 const SummaryState& summary_state,
1790 DeferredLogger& deferred_logger,
1791 double alq_value)
const
1794 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1800 std::vector<double> rates(3);
1801 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1805 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1806 computeBhpAtThpLimitProdWithAlq(frates,
1808 maxPerfPress(ebos_simulator),
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 this->MultisegmentWellGeneric<Scalar>::
1826 computeBhpAtThpLimitProdWithAlq(fratesIter,
1828 maxPerfPress(ebos_simulator),
1838 template<
typename TypeTag>
1839 std::optional<double>
1840 MultisegmentWell<TypeTag>::
1841 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
1842 const SummaryState& summary_state,
1843 DeferredLogger& deferred_logger)
const
1846 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1852 std::vector<double> rates(3);
1853 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1857 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1858 computeBhpAtThpLimitInj(frates,
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 this->MultisegmentWellGeneric<Scalar>::
1876 computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
1883 template<
typename TypeTag>
1885 MultisegmentWell<TypeTag>::
1886 maxPerfPress(
const Simulator& ebos_simulator)
const
1888 double max_pressure = 0.0;
1889 const int nseg = this->numberOfSegments();
1890 for (
int seg = 0; seg < nseg; ++seg) {
1891 for (
const int perf : this->segment_perforations_[seg]) {
1892 const int cell_idx = this->well_cells_[perf];
1893 const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1894 const auto& fs = int_quants.fluidState();
1895 double pressure_cell = this->getPerfCellPressure(fs).value();
1896 max_pressure = std::max(max_pressure, pressure_cell);
1899 return max_pressure;
1906 template<
typename TypeTag>
1913 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
1914 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1915 const int nseg = this->numberOfSegments();
1916 for (
int seg = 0; seg < nseg; ++seg) {
1918 const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
1919 for (
const int perf : this->segment_perforations_[seg]) {
1920 const int cell_idx = this->well_cells_[perf];
1921 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1922 std::vector<Scalar> mob(this->num_components_, 0.0);
1923 getMobilityScalar(ebosSimulator, perf, mob);
1924 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1925 const double Tw = this->well_index_[perf] * trans_mult;
1926 std::vector<Scalar> cq_s(this->num_components_, 0.0);
1927 computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger);
1928 for (
int comp = 0; comp < this->num_components_; ++comp) {
1929 well_q_s[comp] += cq_s[comp];
1940 template<
typename TypeTag>
1944 const std::function<
double(
const double)>& connPICalc,
1945 const std::vector<Scalar>& mobility,
1946 double* connPI)
const
1949 const int np = this->number_of_phases_;
1950 for (
int p = 0; p < np; ++p) {
1953 const auto connMob =
1954 mobility[ this->flowPhaseToEbosCompIdx(p) ]
1955 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
1957 connPI[p] = connPICalc(connMob);
1960 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1961 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1963 const auto io = pu.phase_pos[Oil];
1964 const auto ig = pu.phase_pos[Gas];
1966 const auto vapoil = connPI[ig] * fs.Rv().value();
1967 const auto disgas = connPI[io] * fs.Rs().value();
1969 connPI[io] += vapoil;
1970 connPI[ig] += disgas;
1978 template<
typename TypeTag>
1980 MultisegmentWell<TypeTag>::
1981 computeConnLevelInjInd(
const typename MultisegmentWell<TypeTag>::FluidState& fs,
1982 const Phase preferred_phase,
1983 const std::function<
double(
const double)>& connIICalc,
1984 const std::vector<Scalar>& mobility,
1986 DeferredLogger& deferred_logger)
const
1992 if (preferred_phase == Phase::GAS) {
1993 phase_pos = pu.phase_pos[Gas];
1995 else if (preferred_phase == Phase::OIL) {
1996 phase_pos = pu.phase_pos[Oil];
1998 else if (preferred_phase == Phase::WATER) {
1999 phase_pos = pu.phase_pos[Water];
2002 OPM_DEFLOG_THROW(NotImplemented,
2003 "Unsupported Injector Type ("
2004 <<
static_cast<int>(preferred_phase)
2005 <<
") for well " << this->name()
2006 <<
" during connection I.I. calculation",
2010 const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2011 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:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33