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))
53 , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
56 if constexpr (has_solvent) {
57 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
60 if constexpr (has_polymer) {
61 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
64 if constexpr (Base::has_energy) {
65 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
68 if constexpr (Base::has_foam) {
69 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
72 if constexpr (Base::has_brine) {
73 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
76 if constexpr (Base::has_watVapor) {
77 OPM_THROW(std::runtime_error,
"water evaporation is not supported by multisegment well yet");
80 if(this->rsRvInj() > 0) {
81 OPM_THROW(std::runtime_error,
"dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
82 <<
" \n See (WCONINJE item 10 / WCONHIST item 8)");
90 template <
typename TypeTag>
92 MultisegmentWell<TypeTag>::
93 init(
const PhaseUsage* phase_usage_arg,
94 const std::vector<double>& depth_arg,
95 const double gravity_arg,
97 const std::vector< Scalar >& B_avg,
98 const bool changed_to_open_this_step)
100 Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
112 this->initMatrixAndVectors(num_cells);
115 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
116 const int cell_idx = this->well_cells_[perf];
117 this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
125 template <
typename TypeTag>
127 MultisegmentWell<TypeTag>::
128 initPrimaryVariablesEvaluation()
const
130 this->MSWEval::initPrimaryVariablesEvaluation();
137 template <
typename TypeTag>
139 MultisegmentWell<TypeTag>::
140 updatePrimaryVariables(
const WellState& well_state, DeferredLogger& )
const
142 this->MSWEval::updatePrimaryVariables(well_state);
150 template <
typename TypeTag>
158 Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
161 this->scaleSegmentRatesWithWellRates(well_state);
162 this->scaleSegmentPressuresWithBhp(well_state);
169 template <
typename TypeTag>
173 const std::vector<double>& B_avg,
175 const bool relax_tolerance)
const
177 return this->MSWEval::getWellConvergence(well_state,
180 this->param_.max_residual_allowed_,
181 this->param_.tolerance_wells_,
182 this->param_.relaxed_tolerance_flow_well_,
183 this->param_.tolerance_pressure_ms_wells_,
184 this->param_.relaxed_tolerance_pressure_ms_well_,
192 template <
typename TypeTag>
195 apply(
const BVector& x, BVector& Ax)
const
197 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
199 if ( this->param_.matrix_add_well_contributions_ )
204 BVectorWell Bx(this->duneB_.N());
206 this->duneB_.mv(x, Bx);
209 const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx);
212 this->duneC_.mmtv(invDBx,Ax);
219 template <
typename TypeTag>
222 apply(BVector& r)
const
224 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
227 const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
229 this->duneC_.mmtv(invDrw, r);
234 template <
typename TypeTag>
241 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
244 this->recoverSolutionWell(x, xw);
245 updateWellState(xw, well_state, deferred_logger);
252 template <
typename TypeTag>
257 std::vector<double>& well_potentials,
260 const int np = this->number_of_phases_;
261 well_potentials.resize(np, 0.0);
264 if (this->wellIsStopped()) {
267 this->operability_status_.has_negative_potentials =
false;
270 bool thp_controlled_well =
false;
271 bool bhp_controlled_well =
false;
272 const auto& ws = well_state.well(this->index_of_well_);
273 if (this->isInjector()) {
274 const Well::InjectorCMode& current = ws.injection_cmode;
275 if (current == Well::InjectorCMode::THP) {
276 thp_controlled_well =
true;
278 if (current == Well::InjectorCMode::BHP) {
279 bhp_controlled_well =
true;
282 const Well::ProducerCMode& current = ws.production_cmode;
283 if (current == Well::ProducerCMode::THP) {
284 thp_controlled_well =
true;
286 if (current == Well::ProducerCMode::BHP) {
287 bhp_controlled_well =
true;
290 if (!this->changed_to_open_this_step_ && (thp_controlled_well || bhp_controlled_well)) {
292 double total_rate = 0.0;
293 const double sign = this->isInjector() ? 1.0:-1.0;
294 for (
int phase = 0; phase < np; ++phase){
295 total_rate += sign * ws.surface_rates[phase];
300 if (total_rate > 0) {
301 for (
int phase = 0; phase < np; ++phase){
302 well_potentials[phase] = sign * ws.surface_rates[phase];
308 debug_cost_counter_ = 0;
310 const auto& summaryState = ebosSimulator.vanguard().summaryState();
311 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
312 computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
314 well_potentials = computeWellPotentialWithTHP(
315 well_state, ebosSimulator, deferred_logger);
317 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
318 + this->name() +
": " + std::to_string(debug_cost_counter_));
320 const double sign = this->isInjector() ? 1.0:-1.0;
321 double total_potential = 0.0;
322 for (
int phase = 0; phase < np; ++phase){
323 well_potentials[phase] *= sign;
324 total_potential += well_potentials[phase];
326 if (total_potential < 0.0 && this->param_.check_well_operability_) {
328 this->operability_status_.has_negative_potentials =
true;
329 const std::string msg = std::string(
"well ") + this->name() + std::string(
": has non negative potentials is not operable");
330 deferred_logger.warning(
"NEGATIVE_POTENTIALS_INOPERABLE", msg);
337 template<
typename TypeTag>
341 std::vector<double>& well_flux,
344 if (this->well_ecl_.isInjector()) {
345 const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
346 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
348 const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
349 computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
353 template<
typename TypeTag>
355 MultisegmentWell<TypeTag>::
356 computeWellRatesWithBhp(
const Simulator& ebosSimulator,
358 std::vector<double>& well_flux,
359 DeferredLogger& deferred_logger)
const
362 const int np = this->number_of_phases_;
364 well_flux.resize(np, 0.0);
365 const bool allow_cf = this->getAllowCrossFlow();
366 const int nseg = this->numberOfSegments();
367 const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
368 const auto& ws = well_state.well(this->indexOfWell());
369 auto segments_copy = ws.segments;
370 segments_copy.scale_pressure(bhp);
371 const auto& segment_pressure = segments_copy.pressure;
372 for (
int seg = 0; seg < nseg; ++seg) {
373 for (
const int perf : this->segment_perforations_[seg]) {
374 const int cell_idx = this->well_cells_[perf];
375 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
377 std::vector<Scalar> mob(this->num_components_, 0.);
378 getMobilityScalar(ebosSimulator, perf, mob);
379 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
380 const double Tw = this->well_index_[perf] * trans_mult;
382 const Scalar seg_pressure = segment_pressure[seg];
383 std::vector<Scalar> cq_s(this->num_components_, 0.);
384 computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
385 allow_cf, cq_s, deferred_logger);
387 for(
int p = 0; p < np; ++p) {
388 well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
392 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
396 template<
typename TypeTag>
398 MultisegmentWell<TypeTag>::
399 computeWellRatesWithBhpIterations(
const Simulator& ebosSimulator,
401 std::vector<double>& well_flux,
402 DeferredLogger& deferred_logger)
const
406 MultisegmentWell<TypeTag> well_copy(*
this);
407 well_copy.debug_cost_counter_ = 0;
410 WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
411 const auto& group_state = ebosSimulator.problem().wellModel().groupState();
412 auto& ws = well_state_copy.well(this->index_of_well_);
415 const auto& summary_state = ebosSimulator.vanguard().summaryState();
416 auto inj_controls = well_copy.well_ecl_.isInjector()
417 ? well_copy.well_ecl_.injectionControls(summary_state)
418 : Well::InjectionControls(0);
419 auto prod_controls = well_copy.well_ecl_.isProducer()
420 ? well_copy.well_ecl_.productionControls(summary_state) :
421 Well::ProductionControls(0);
424 if (well_copy.well_ecl_.isInjector()) {
425 inj_controls.bhp_limit = bhp;
426 ws.injection_cmode = Well::InjectorCMode::BHP;
428 prod_controls.bhp_limit = bhp;
429 ws.production_cmode = Well::ProducerCMode::BHP;
432 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
435 const int np = this->number_of_phases_;
437 for (
int phase = 0; phase < np; ++phase){
438 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
441 const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
442 for (
int phase = 0; phase < np; ++phase) {
443 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
446 well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
448 well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
449 const double dt = ebosSimulator.timeStepSize();
451 well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
456 well_flux.resize(np, 0.0);
457 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
458 const EvalWell rate = well_copy.getQs(compIdx);
459 well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
461 debug_cost_counter_ += well_copy.debug_cost_counter_;
466 template<
typename TypeTag>
468 MultisegmentWell<TypeTag>::
469 computeWellPotentialWithTHP(
470 const WellState& well_state,
471 const Simulator& ebos_simulator,
472 DeferredLogger& deferred_logger)
const
474 std::vector<double> potentials(this->number_of_phases_, 0.0);
475 const auto& summary_state = ebos_simulator.vanguard().summaryState();
477 const auto& well = this->well_ecl_;
478 if (well.isInjector()){
479 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
480 if (bhp_at_thp_limit) {
481 const auto& controls = well.injectionControls(summary_state);
482 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
483 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
484 deferred_logger.debug(
"Converged thp based potential calculation for well "
485 + this->name() +
", at bhp = " + std::to_string(bhp));
487 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
488 "Failed in getting converged thp based potential calculation for well "
489 + this->name() +
". Instead the bhp based value is used");
490 const auto& controls = well.injectionControls(summary_state);
491 const double bhp = controls.bhp_limit;
492 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
495 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
496 well_state, ebos_simulator, summary_state, deferred_logger);
497 if (bhp_at_thp_limit) {
498 const auto& controls = well.productionControls(summary_state);
499 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
500 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
501 deferred_logger.debug(
"Converged thp based potential calculation for well "
502 + this->name() +
", at bhp = " + std::to_string(bhp));
504 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
505 "Failed in getting converged thp based potential calculation for well "
506 + this->name() +
". Instead the bhp based value is used");
507 const auto& controls = well.productionControls(summary_state);
508 const double bhp = controls.bhp_limit;
509 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
518 template <
typename TypeTag>
520 MultisegmentWell<TypeTag>::
521 solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
523 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
527 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
529 updateWellState(dx_well, well_state, deferred_logger);
536 template <
typename TypeTag>
538 MultisegmentWell<TypeTag>::
539 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
541 for (
int perf = 0; perf < this->number_of_perforations_; ++perf) {
543 std::vector<double> kr(this->number_of_phases_, 0.0);
544 std::vector<double> density(this->number_of_phases_, 0.0);
546 const int cell_idx = this->well_cells_[perf];
547 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
548 const auto& fs = intQuants.fluidState();
553 if (pu.phase_used[Water]) {
554 const int water_pos = pu.phase_pos[Water];
555 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
556 sum_kr += kr[water_pos];
557 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
560 if (pu.phase_used[Oil]) {
561 const int oil_pos = pu.phase_pos[Oil];
562 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
563 sum_kr += kr[oil_pos];
564 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
567 if (pu.phase_used[Gas]) {
568 const int gas_pos = pu.phase_pos[Gas];
569 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
570 sum_kr += kr[gas_pos];
571 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
574 assert(sum_kr != 0.);
577 double average_density = 0.;
578 for (
int p = 0; p < this->number_of_phases_; ++p) {
579 average_density += kr[p] * density[p];
581 average_density /= sum_kr;
583 this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
591 template <
typename TypeTag>
593 MultisegmentWell<TypeTag>::
594 computeInitialSegmentFluids(
const Simulator& ebos_simulator)
596 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
598 const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
599 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
600 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
609 template <
typename TypeTag>
611 MultisegmentWell<TypeTag>::
612 updateWellState(
const BVectorWell& dwells,
613 WellState& well_state,
614 DeferredLogger& deferred_logger,
615 const double relaxation_factor)
const
617 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
619 const double dFLimit = this->param_.dwell_fraction_max_;
620 const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
621 this->MSWEval::updatePrimaryVariablesNewton(dwells,
624 max_pressure_change);
626 this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
627 Base::calculateReservoirRates(well_state.well(this->index_of_well_));
634 template <
typename TypeTag>
636 MultisegmentWell<TypeTag>::
637 calculateExplicitQuantities(
const Simulator& ebosSimulator,
638 const WellState& well_state,
639 DeferredLogger& deferred_logger)
641 updatePrimaryVariables(well_state, deferred_logger);
642 initPrimaryVariablesEvaluation();
643 computePerfCellPressDiffs(ebosSimulator);
644 computeInitialSegmentFluids(ebosSimulator);
651 template<
typename TypeTag>
653 MultisegmentWell<TypeTag>::
654 updateProductivityIndex(
const Simulator& ebosSimulator,
655 const WellProdIndexCalculator& wellPICalc,
656 WellState& well_state,
657 DeferredLogger& deferred_logger)
const
659 auto fluidState = [&ebosSimulator,
this](
const int perf)
661 const auto cell_idx = this->well_cells_[perf];
662 return ebosSimulator.model()
663 .cachedIntensiveQuantities(cell_idx, 0)->fluidState();
666 const int np = this->number_of_phases_;
667 auto setToZero = [np](
double* x) ->
void
669 std::fill_n(x, np, 0.0);
672 auto addVector = [np](
const double* src,
double* dest) ->
void
674 std::transform(src, src + np, dest, dest, std::plus<>{});
677 auto& ws = well_state.well(this->index_of_well_);
678 auto& perf_data = ws.perf_data;
679 auto* connPI = perf_data.prod_index.data();
680 auto* wellPI = ws.productivity_index.data();
684 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
685 auto subsetPerfID = 0;
687 for (
const auto& perf : *this->perf_data_){
688 auto allPerfID = perf.ecl_index;
690 auto connPICalc = [&wellPICalc, allPerfID](
const double mobility) ->
double
692 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
695 std::vector<Scalar> mob(this->num_components_, 0.0);
696 getMobilityScalar(ebosSimulator,
static_cast<int>(subsetPerfID), mob);
698 const auto& fs = fluidState(subsetPerfID);
701 if (this->isInjector()) {
702 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
703 mob, connPI, deferred_logger);
706 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
709 addVector(connPI, wellPI);
715 assert (
static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
716 "Internal logic error in processing connections for PI/II");
723 template<
typename TypeTag>
725 MultisegmentWell<TypeTag>::
726 addWellContributions(SparseMatrixAdapter& jacobian)
const
728 const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(this->duneD_, this->duneDSolver_);
738 for (
size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
739 for (
auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
740 const auto row_index = colC.index();
741 for (
size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
742 for (
auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
743 const auto col_index = colB.index();
744 OffDiagMatrixBlockWellType tmp1;
745 detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
746 typename SparseMatrixAdapter::MatrixBlock tmp2;
747 detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
748 jacobian.addToBlock(row_index, col_index, tmp2);
756 template<
typename TypeTag>
758 MultisegmentWell<TypeTag>::
759 addWellPressureEquations(PressureMatrix& jacobian,
760 const BVector& weights,
761 const int pressureVarIndex,
763 const WellState& well_state)
const
768 const auto seg_pressure_var_ind = this->SPres;
769 const int welldof_ind = this->duneC_.M() + this->index_of_well_;
770 if(not(this->isPressureControlled(well_state))){
771 for (
size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
772 for (
auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
773 const auto row_index = colC.index();
774 const auto& bw = weights[row_index];
777 for(
size_t i = 0; i< bw.size(); ++i){
778 matel += bw[i]*(*colC)[seg_pressure_var_ind][i];
780 jacobian[row_index][welldof_ind] += matel;
786 if(not(this->isPressureControlled(well_state))){
787 auto well_weight = weights[0];
790 for (
size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
791 for (
auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
792 const auto col_index = colB.index();
793 const auto& bw = weights[col_index];
799 well_weight /= num_perfs;
803 double diag_ell = 0.0;
804 for (
size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
805 const auto& bw = well_weight;
806 for (
auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
807 const auto col_index = colB.index();
809 for(
size_t i = 0; i< bw.size(); ++i){
810 matel += bw[i] *(*colB)[i][pressureVarIndex];
812 jacobian[welldof_ind][col_index] += matel;
817 assert(diag_ell > 0.0);
818 jacobian[welldof_ind][welldof_ind] = diag_ell;
820 jacobian[welldof_ind][welldof_ind] = 1.0;
825 template<
typename TypeTag>
826 template<
class Value>
828 MultisegmentWell<TypeTag>::
829 computePerfRate(
const Value& pressure_cell,
832 const std::vector<Value>& b_perfcells,
833 const std::vector<Value>& mob_perfcells,
836 const Value& segment_pressure,
837 const Value& segment_density,
838 const bool& allow_cf,
839 const std::vector<Value>& cmix_s,
840 std::vector<Value>& cq_s,
842 double& perf_dis_gas_rate,
843 double& perf_vap_oil_rate,
844 DeferredLogger& deferred_logger)
const
847 const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
849 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
851 perf_press = pressure_cell - cell_perf_press_diff;
854 const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
857 if ( drawdown > 0.0) {
859 if (!allow_cf && this->isInjector()) {
864 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
865 const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
866 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
869 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
870 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
871 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
872 const Value cq_s_oil = cq_s[oilCompIdx];
873 const Value cq_s_gas = cq_s[gasCompIdx];
874 cq_s[gasCompIdx] += rs * cq_s_oil;
875 cq_s[oilCompIdx] += rv * cq_s_gas;
879 if (!allow_cf && this->isProducer()) {
884 Value total_mob = mob_perfcells[0];
885 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
886 total_mob += mob_perfcells[comp_idx];
890 const Value cqt_i = - Tw * (total_mob * drawdown);
893 Value volume_ratio = 0.0;
894 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
895 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
896 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
899 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
900 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
901 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
906 const Value d = 1.0 - rv * rs;
908 if (getValue(d) == 0.0) {
909 OPM_DEFLOG_THROW(NumericalIssue,
"Zero d value obtained for well " << this->name()
910 <<
" during flux calculation"
911 <<
" with rs " << rs <<
" and rv " << rv, deferred_logger);
914 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
915 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
917 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
918 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
920 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
921 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
922 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
924 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
925 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
926 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
930 Value cqt_is = cqt_i / volume_ratio;
931 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
932 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
937 if (this->isProducer()) {
938 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
939 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
940 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
949 const double d = 1.0 - getValue(rv) * getValue(rs);
952 perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
955 perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
960 template <
typename TypeTag>
962 MultisegmentWell<TypeTag>::
963 computePerfRateEval(
const IntensiveQuantities& int_quants,
964 const std::vector<EvalWell>& mob_perfcells,
968 const EvalWell& segment_pressure,
969 const bool& allow_cf,
970 std::vector<EvalWell>& cq_s,
971 EvalWell& perf_press,
972 double& perf_dis_gas_rate,
973 double& perf_vap_oil_rate,
974 DeferredLogger& deferred_logger)
const
977 const auto& fs = int_quants.fluidState();
979 const EvalWell pressure_cell = this->extendEval(this->getPerfCellPressure(fs));
980 const EvalWell rs = this->extendEval(fs.Rs());
981 const EvalWell rv = this->extendEval(fs.Rv());
984 std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
986 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
987 if (!FluidSystem::phaseIsActive(phaseIdx)) {
991 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
992 b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
995 std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
996 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
997 cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
1000 this->computePerfRate(pressure_cell,
1008 this->segment_densities_[seg],
1020 template <
typename TypeTag>
1022 MultisegmentWell<TypeTag>::
1023 computePerfRateScalar(
const IntensiveQuantities& int_quants,
1024 const std::vector<Scalar>& mob_perfcells,
1028 const Scalar& segment_pressure,
1029 const bool& allow_cf,
1030 std::vector<Scalar>& cq_s,
1031 DeferredLogger& deferred_logger)
const
1034 const auto& fs = int_quants.fluidState();
1036 const Scalar pressure_cell = getValue(this->getPerfCellPressure(fs));
1037 const Scalar rs = getValue(fs.Rs());
1038 const Scalar rv = getValue(fs.Rv());
1041 std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
1043 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1044 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1048 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1049 b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
1052 std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
1053 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
1054 cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
1057 Scalar perf_dis_gas_rate = 0.0;
1058 Scalar perf_vap_oil_rate = 0.0;
1059 Scalar perf_press = 0.0;
1061 this->computePerfRate(pressure_cell,
1069 getValue(this->segment_densities_[seg]),
1079 template <
typename TypeTag>
1081 MultisegmentWell<TypeTag>::
1082 computeSegmentFluidProperties(
const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
1091 EvalWell temperature;
1092 EvalWell saltConcentration;
1098 int pvt_region_index;
1101 const int cell_idx = this->well_cells_[0];
1102 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1103 const auto& fs = intQuants.fluidState();
1104 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1105 saltConcentration = this->extendEval(fs.saltConcentration());
1106 pvt_region_index = fs.pvtRegionIndex();
1109 this->MSWEval::computeSegmentFluidProperties(temperature,
1119 template <
typename TypeTag>
1121 MultisegmentWell<TypeTag>::
1122 getMobilityEval(
const Simulator& ebosSimulator,
1124 std::vector<EvalWell>& mob)
const
1127 const int cell_idx = this->well_cells_[perf];
1128 assert (
int(mob.size()) == this->num_components_);
1129 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1130 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1134 const int satid = this->saturation_table_number_[perf] - 1;
1135 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1136 if( satid == satid_elem ) {
1138 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1139 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1143 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1144 mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
1151 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1152 std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
1153 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1156 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1159 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1160 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1164 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1165 mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1171 template <
typename TypeTag>
1173 MultisegmentWell<TypeTag>::
1174 getMobilityScalar(
const Simulator& ebosSimulator,
1176 std::vector<Scalar>& mob)
const
1179 const int cell_idx = this->well_cells_[perf];
1180 assert (
int(mob.size()) == this->num_components_);
1181 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1182 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1186 const int satid = this->saturation_table_number_[perf] - 1;
1187 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1188 if( satid == satid_elem ) {
1190 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1191 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1195 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1196 mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
1203 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1204 std::array<Scalar,3> relativePerms = { 0.0, 0.0, 0.0 };
1205 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1208 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1211 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1212 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1216 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1217 mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
1225 template<
typename TypeTag>
1227 MultisegmentWell<TypeTag>::
1228 getRefDensity()
const
1230 return this->segment_densities_[0].value();
1233 template<
typename TypeTag>
1235 MultisegmentWell<TypeTag>::
1236 checkOperabilityUnderBHPLimit(
const WellState& ,
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1238 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1239 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1242 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1243 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1246 double total_ipr_mass_rate = 0.0;
1247 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1249 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1253 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1254 const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1256 const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1257 total_ipr_mass_rate += ipr_rate * rho;
1259 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1260 this->operability_status_.operable_under_only_bhp_limit =
false;
1264 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1268 std::vector<double> well_rates_bhp_limit;
1269 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1271 const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
1272 const double thp_limit = this->getTHPConstraint(summaryState);
1273 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1274 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1285 this->operability_status_.operable_under_only_bhp_limit =
true;
1286 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1292 template<
typename TypeTag>
1294 MultisegmentWell<TypeTag>::
1295 updateIPR(
const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
const
1300 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1301 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1303 const int nseg = this->numberOfSegments();
1304 const double ref_depth = this->ref_depth_;
1305 std::vector<double> seg_dp(nseg, 0.0);
1306 for (
int seg = 0; seg < nseg; ++seg) {
1308 const double segment_depth = this->segmentSet()[seg].depth();
1309 const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
1310 const double segment_depth_outlet = seg == 0? ref_depth : this->segmentSet()[outlet_segment_index].depth();
1311 double dp = wellhelpers::computeHydrostaticCorrection(segment_depth_outlet, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
1315 dp += seg_dp[outlet_segment_index];
1318 for (
const int perf : this->segment_perforations_[seg]) {
1319 std::vector<Scalar> mob(this->num_components_, 0.0);
1322 getMobilityScalar(ebos_simulator, perf, mob);
1324 const int cell_idx = this->well_cells_[perf];
1325 const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1326 const auto& fs = int_quantities.fluidState();
1329 const double perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
1331 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1332 const double pressure_cell = this->getPerfCellPressure(fs).value();
1335 std::vector<double> b_perf(this->num_components_);
1336 for (
size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1337 if (!FluidSystem::phaseIsActive(phase)) {
1340 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1341 b_perf[comp_idx] = fs.invB(phase).value();
1345 const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1346 const double pressure_diff = pressure_cell - h_perf;
1349 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1350 deferred_logger.debug(
"CROSSFLOW_IPR",
1351 "cross flow found when updateIPR for well " + this->name());
1355 const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1357 std::vector<double> ipr_a_perf(this->ipr_a_.size());
1358 std::vector<double> ipr_b_perf(this->ipr_b_.size());
1359 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1360 const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1361 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1362 ipr_b_perf[comp_idx] += tw_mob;
1366 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1367 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1368 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1369 const double rs = (fs.Rs()).value();
1370 const double rv = (fs.Rv()).value();
1372 const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1373 const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1375 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1376 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1378 const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1379 const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1381 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1382 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1385 for (
size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1386 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1387 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1393 template<
typename TypeTag>
1395 MultisegmentWell<TypeTag>::
1396 checkOperabilityUnderTHPLimit(
1397 const Simulator& ebos_simulator,
1398 const WellState& well_state,
1399 DeferredLogger& deferred_logger)
1401 const auto& summaryState = ebos_simulator.vanguard().summaryState();
1402 const auto obtain_bhp = this->isProducer()
1403 ? computeBhpAtThpLimitProd(
1404 well_state, ebos_simulator, summaryState, deferred_logger)
1405 : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1408 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1410 const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1411 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1413 const double thp_limit = this->getTHPConstraint(summaryState);
1414 if (this->isProducer() && *obtain_bhp < thp_limit) {
1415 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1416 +
" bars is SMALLER than thp limit "
1417 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1418 +
" bars as a producer for well " + this->name();
1419 deferred_logger.debug(msg);
1421 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1422 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1423 +
" bars is LARGER than thp limit "
1424 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1425 +
" bars as a injector for well " + this->name();
1426 deferred_logger.debug(msg);
1431 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1432 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1433 if (!this->wellIsStopped()) {
1434 const double thp_limit = this->getTHPConstraint(summaryState);
1435 deferred_logger.debug(
" could not find bhp value at thp limit "
1436 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1437 +
" bar for well " + this->name() +
", the well might need to be closed ");
1446 template<
typename TypeTag>
1448 MultisegmentWell<TypeTag>::
1449 iterateWellEqWithControl(
const Simulator& ebosSimulator,
1451 const Well::InjectionControls& inj_controls,
1452 const Well::ProductionControls& prod_controls,
1453 WellState& well_state,
1454 const GroupState& group_state,
1455 DeferredLogger& deferred_logger)
1457 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1459 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1463 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1468 std::vector<std::vector<Scalar> > residual_history;
1469 std::vector<double> measure_history;
1472 double relaxation_factor = 1.;
1473 const double min_relaxation_factor = 0.6;
1474 bool converged =
false;
1475 int stagnate_count = 0;
1476 bool relax_convergence =
false;
1477 this->regularize_ =
false;
1478 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1480 assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1482 const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
1484 if (it > this->param_.strict_inner_iter_wells_) {
1485 relax_convergence =
true;
1486 this->regularize_ =
true;
1489 const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
1490 if (report.converged()) {
1497 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1501 residual_history.push_back(residuals);
1502 measure_history.push_back(this->getResidualMeasureValue(well_state,
1503 residual_history[it],
1504 this->param_.tolerance_wells_,
1505 this->param_.tolerance_pressure_ms_wells_,
1510 bool is_oscillate =
false;
1511 bool is_stagnate =
false;
1513 this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
1517 if (is_oscillate || is_stagnate) {
1519 std::ostringstream sstr;
1520 if (relaxation_factor == min_relaxation_factor) {
1523 if (stagnate_count == 6) {
1524 sstr <<
" well " << this->name() <<
" observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1525 const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger,
true);
1526 if (reportStag.converged()) {
1528 sstr <<
" well " << this->name() <<
" manages to get converged with relaxed tolerances in " << it <<
" inner iterations";
1529 deferred_logger.debug(sstr.str());
1536 const double reduction_mutliplier = 0.9;
1537 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1541 sstr <<
" well " << this->name() <<
" observes stagnation in inner iteration " << it <<
"\n";
1545 sstr <<
" well " << this->name() <<
" observes oscillation in inner iteration " << it <<
"\n";
1547 sstr <<
" relaxation_factor is " << relaxation_factor <<
" now\n";
1549 this->regularize_ =
true;
1550 deferred_logger.debug(sstr.str());
1552 updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
1553 initPrimaryVariablesEvaluation();
1558 std::ostringstream sstr;
1559 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1560 if (relax_convergence)
1561 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1562 deferred_logger.debug(sstr.str());
1564 std::ostringstream sstr;
1565 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1566#define EXTRA_DEBUG_MSW 0
1568 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1569 for (
int i = 0; i < it; ++i) {
1570 const auto& residual = residual_history[i];
1571 sstr <<
" residual at " << i <<
"th iteration ";
1572 for (
const auto& res : residual) {
1575 sstr <<
" " << measure_history[i] <<
" \n";
1578 deferred_logger.debug(sstr.str());
1588 template<
typename TypeTag>
1590 MultisegmentWell<TypeTag>::
1591 assembleWellEqWithoutIteration(
const Simulator& ebosSimulator,
1593 const Well::InjectionControls& inj_controls,
1594 const Well::ProductionControls& prod_controls,
1595 WellState& well_state,
1596 const GroupState& group_state,
1597 DeferredLogger& deferred_logger)
1600 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1603 this->updateUpwindingSegments();
1606 computeSegmentFluidProperties(ebosSimulator, deferred_logger);
1613 this->resWell_ = 0.0;
1615 this->duneDSolver_.reset();
1617 auto& ws = well_state.well(this->index_of_well_);
1618 ws.dissolved_gas_rate = 0;
1619 ws.vaporized_oil_rate = 0;
1620 ws.vaporized_wat_rate = 0;
1627 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1629 const int nseg = this->numberOfSegments();
1631 for (
int seg = 0; seg < nseg; ++seg) {
1635 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1640 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1642 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1643 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
1644 - segment_fluid_initial_[seg][comp_idx]) / dt;
1646 this->resWell_[seg][comp_idx] += accumulation_term.value();
1647 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1648 this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
1654 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1655 const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
1657 const int seg_upwind = this->upwinding_segments_[seg];
1660 this->resWell_[seg][comp_idx] -= segment_rate.value();
1661 this->duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
1662 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1663 this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
1665 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1666 this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
1674 for (
const int inlet : this->segment_inlets_[seg]) {
1675 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1676 const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
1678 const int inlet_upwind = this->upwinding_segments_[inlet];
1681 this->resWell_[seg][comp_idx] += inlet_rate.value();
1682 this->duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
1683 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1684 this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
1686 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1687 this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
1695 const EvalWell seg_pressure = this->getSegmentPressure(seg);
1696 auto& perf_data = ws.perf_data;
1697 auto& perf_rates = perf_data.phase_rates;
1698 auto& perf_press_state = perf_data.pressure;
1699 for (
const int perf : this->segment_perforations_[seg]) {
1700 const int cell_idx = this->well_cells_[perf];
1701 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1702 std::vector<EvalWell> mob(this->num_components_, 0.0);
1703 getMobilityEval(ebosSimulator, perf, mob);
1704 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1705 const double Tw = this->well_index_[perf] * trans_mult;
1706 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1707 EvalWell perf_press;
1708 double perf_dis_gas_rate = 0.;
1709 double perf_vap_oil_rate = 0.;
1710 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);
1713 if (this->isProducer()) {
1714 ws.dissolved_gas_rate += perf_dis_gas_rate;
1715 ws.vaporized_oil_rate += perf_vap_oil_rate;
1719 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1720 perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1722 perf_press_state[perf] = perf_press.value();
1724 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1726 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1728 this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1731 this->resWell_[seg][comp_idx] += cq_s_effective.value();
1734 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1737 this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq);
1740 this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
1743 for (
int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
1745 this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
1752 const auto& summaryState = ebosSimulator.vanguard().summaryState();
1753 const Schedule& schedule = ebosSimulator.vanguard().schedule();
1754 this->assembleControlEq(well_state,
1763 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1764 this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
1772 template<
typename TypeTag>
1774 MultisegmentWell<TypeTag>::
1775 openCrossFlowAvoidSingularity(
const Simulator& ebos_simulator)
const
1777 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1781 template<
typename TypeTag>
1783 MultisegmentWell<TypeTag>::
1784 allDrawDownWrongDirection(
const Simulator& ebos_simulator)
const
1786 bool all_drawdown_wrong_direction =
true;
1787 const int nseg = this->numberOfSegments();
1789 for (
int seg = 0; seg < nseg; ++seg) {
1790 const EvalWell segment_pressure = this->getSegmentPressure(seg);
1791 for (
const int perf : this->segment_perforations_[seg]) {
1793 const int cell_idx = this->well_cells_[perf];
1794 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1795 const auto& fs = intQuants.fluidState();
1798 const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
1800 const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1802 const double pressure_cell = this->getPerfCellPressure(fs).value();
1803 const double perf_press = pressure_cell - cell_perf_press_diff;
1806 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1811 if ( (drawdown < 0. && this->isInjector()) ||
1812 (drawdown > 0. && this->isProducer()) ) {
1813 all_drawdown_wrong_direction =
false;
1819 return all_drawdown_wrong_direction;
1825 template<
typename TypeTag>
1827 MultisegmentWell<TypeTag>::
1828 updateWaterThroughput(
const double , WellState& )
const
1836 template<
typename TypeTag>
1837 typename MultisegmentWell<TypeTag>::EvalWell
1838 MultisegmentWell<TypeTag>::
1839 getSegmentSurfaceVolume(
const Simulator& ebos_simulator,
const int seg_idx)
const
1841 EvalWell temperature;
1842 EvalWell saltConcentration;
1843 int pvt_region_index;
1847 const int cell_idx = this->well_cells_[0];
1848 const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1849 const auto& fs = intQuants.fluidState();
1850 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1851 saltConcentration = this->extendEval(fs.saltConcentration());
1852 pvt_region_index = fs.pvtRegionIndex();
1855 return this->MSWEval::getSegmentSurfaceVolume(temperature,
1862 template<
typename TypeTag>
1863 std::optional<double>
1864 MultisegmentWell<TypeTag>::
1865 computeBhpAtThpLimitProd(
const WellState& well_state,
1866 const Simulator& ebos_simulator,
1867 const SummaryState& summary_state,
1868 DeferredLogger& deferred_logger)
const
1870 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
1874 this->getALQ(well_state));
1879 template<
typename TypeTag>
1880 std::optional<double>
1881 MultisegmentWell<TypeTag>::
1882 computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
1883 const SummaryState& summary_state,
1884 DeferredLogger& deferred_logger,
1885 double alq_value)
const
1888 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1894 std::vector<double> rates(3);
1895 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1899 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1900 computeBhpAtThpLimitProdWithAlq(frates,
1902 maxPerfPress(ebos_simulator),
1910 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1914 std::vector<double> rates(3);
1915 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1919 return this->MultisegmentWellGeneric<Scalar>::
1920 computeBhpAtThpLimitProdWithAlq(fratesIter,
1922 maxPerfPress(ebos_simulator),
1932 template<
typename TypeTag>
1933 std::optional<double>
1934 MultisegmentWell<TypeTag>::
1935 computeBhpAtThpLimitInj(
const Simulator& ebos_simulator,
1936 const SummaryState& summary_state,
1937 DeferredLogger& deferred_logger)
const
1940 auto frates = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1946 std::vector<double> rates(3);
1947 computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1951 auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1952 computeBhpAtThpLimitInj(frates,
1960 auto fratesIter = [
this, &ebos_simulator, &deferred_logger](
const double bhp) {
1964 std::vector<double> rates(3);
1965 computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1969 return this->MultisegmentWellGeneric<Scalar>::
1970 computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
1977 template<
typename TypeTag>
1979 MultisegmentWell<TypeTag>::
1980 maxPerfPress(
const Simulator& ebos_simulator)
const
1982 double max_pressure = 0.0;
1983 const int nseg = this->numberOfSegments();
1984 for (
int seg = 0; seg < nseg; ++seg) {
1985 for (
const int perf : this->segment_perforations_[seg]) {
1986 const int cell_idx = this->well_cells_[perf];
1987 const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, 0));
1988 const auto& fs = int_quants.fluidState();
1989 double pressure_cell = this->getPerfCellPressure(fs).value();
1990 max_pressure = std::max(max_pressure, pressure_cell);
1993 return max_pressure;
2000 template<
typename TypeTag>
2007 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
2008 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2009 const int nseg = this->numberOfSegments();
2010 for (
int seg = 0; seg < nseg; ++seg) {
2012 const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
2013 for (
const int perf : this->segment_perforations_[seg]) {
2014 const int cell_idx = this->well_cells_[perf];
2015 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
2016 std::vector<Scalar> mob(this->num_components_, 0.0);
2017 getMobilityScalar(ebosSimulator, perf, mob);
2018 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
2019 const double Tw = this->well_index_[perf] * trans_mult;
2020 std::vector<Scalar> cq_s(this->num_components_, 0.0);
2021 computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger);
2022 for (
int comp = 0; comp < this->num_components_; ++comp) {
2023 well_q_s[comp] += cq_s[comp];
2034 template<
typename TypeTag>
2038 const std::function<
double(
const double)>& connPICalc,
2039 const std::vector<Scalar>& mobility,
2040 double* connPI)
const
2043 const int np = this->number_of_phases_;
2044 for (
int p = 0; p < np; ++p) {
2047 const auto connMob =
2048 mobility[ this->flowPhaseToEbosCompIdx(p) ]
2049 * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
2051 connPI[p] = connPICalc(connMob);
2054 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2055 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2057 const auto io = pu.phase_pos[Oil];
2058 const auto ig = pu.phase_pos[Gas];
2060 const auto vapoil = connPI[ig] * fs.Rv().value();
2061 const auto disgas = connPI[io] * fs.Rs().value();
2063 connPI[io] += vapoil;
2064 connPI[ig] += disgas;
2072 template<
typename TypeTag>
2074 MultisegmentWell<TypeTag>::
2075 computeConnLevelInjInd(
const typename MultisegmentWell<TypeTag>::FluidState& fs,
2076 const Phase preferred_phase,
2077 const std::function<
double(
const double)>& connIICalc,
2078 const std::vector<Scalar>& mobility,
2080 DeferredLogger& deferred_logger)
const
2086 if (preferred_phase == Phase::GAS) {
2087 phase_pos = pu.phase_pos[Gas];
2089 else if (preferred_phase == Phase::OIL) {
2090 phase_pos = pu.phase_pos[Oil];
2092 else if (preferred_phase == Phase::WATER) {
2093 phase_pos = pu.phase_pos[Water];
2096 OPM_DEFLOG_THROW(NotImplemented,
2097 "Unsupported Injector Type ("
2098 <<
static_cast<int>(preferred_phase)
2099 <<
") for well " << this->name()
2100 <<
" during connection I.I. calculation",
2104 const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
2105 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:37