22#ifndef OPM_MULTISEGMENTWELL_IMPL_HEADER_INCLUDED
23#define OPM_MULTISEGMENTWELL_IMPL_HEADER_INCLUDED
25#ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
27#include <opm/simulators/wells/MultisegmentWell.hpp>
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/OpmLog/OpmLog.hpp>
33#include <opm/input/eclipse/Schedule/MSW/Segment.hpp>
34#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
35#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
36#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
37#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
39#include <opm/input/eclipse/Units/Units.hpp>
41#include <opm/material/densead/EvaluationFormat.hpp>
43#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
44#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
45#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
51#if COMPILE_GPU_BRIDGE && (HAVE_CUDA || HAVE_OPENCL)
52#include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
59 template <
typename TypeTag>
60 MultisegmentWell<TypeTag>::
61 MultisegmentWell(
const Well& well,
64 const ModelParameters& param,
66 const int pvtRegionIdx,
74 , segment_fluid_initial_(
this->numberOfSegments(), std::
vector<Scalar>(
this->num_components_, 0.0))
77 if constexpr (has_solvent) {
78 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
81 if constexpr (has_polymer) {
82 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
85 if constexpr (Base::has_energy) {
86 OPM_THROW(std::runtime_error,
"energy is not supported by multisegment well yet");
89 if constexpr (Base::has_foam) {
90 OPM_THROW(std::runtime_error,
"foam is not supported by multisegment well yet");
93 if constexpr (Base::has_brine) {
94 OPM_THROW(std::runtime_error,
"brine is not supported by multisegment well yet");
97 if constexpr (Base::has_watVapor) {
98 OPM_THROW(std::runtime_error,
"water evaporation is not supported by multisegment well yet");
101 if constexpr (Base::has_micp) {
102 OPM_THROW(std::runtime_error,
"MICP is not supported by multisegment well yet");
105 if(this->rsRvInj() > 0) {
106 OPM_THROW(std::runtime_error,
107 "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
108 " \n See (WCONINJE item 10 / WCONHIST item 8)");
111 this->thp_update_iterations =
true;
118 template <
typename TypeTag>
120 MultisegmentWell<TypeTag>::
121 init(
const PhaseUsage* phase_usage_arg,
122 const std::vector<Scalar>& depth_arg,
123 const Scalar gravity_arg,
124 const std::vector< Scalar >& B_avg,
125 const bool changed_to_open_this_step)
127 Base::init(phase_usage_arg, depth_arg, gravity_arg, B_avg, changed_to_open_this_step);
139 this->initMatrixAndVectors();
142 for (
int local_perf_index = 0; local_perf_index < this->number_of_local_perforations_; ++local_perf_index) {
144 const int cell_idx = this->well_cells_[local_perf_index];
146 this->cell_perforation_depth_diffs_[local_perf_index] = depth_arg[cell_idx] - this->perf_depth_[this->pw_info_.localToActive(local_perf_index)];
154 template <
typename TypeTag>
156 MultisegmentWell<TypeTag>::
157 updatePrimaryVariables(
const Simulator& simulator,
158 const WellState<Scalar>& well_state,
159 DeferredLogger& deferred_logger)
161 const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
162 this->primary_variables_.update(well_state, stop_or_zero_rate_target);
170 template <
typename TypeTag>
178 Base::updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
181 this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
182 this->segments_.perforations(),
184 this->scaleSegmentPressuresWithBhp(well_state);
191 template <
typename TypeTag>
196 const std::vector<Scalar>& B_avg,
198 const bool relax_tolerance)
const
200 return this->MSWEval::getWellConvergence(well_state,
203 this->param_.max_residual_allowed_,
204 this->param_.tolerance_wells_,
205 this->param_.relaxed_tolerance_flow_well_,
206 this->param_.tolerance_pressure_ms_wells_,
207 this->param_.relaxed_tolerance_pressure_ms_well_,
209 this->wellIsStopped());
217 template <
typename TypeTag>
220 apply(
const BVector& x, BVector& Ax)
const
222 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
226 if (this->param_.matrix_add_well_contributions_) {
231 this->linSys_.
apply(x, Ax);
238 template <
typename TypeTag>
241 apply(BVector& r)
const
243 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
247 this->linSys_.
apply(r);
252 template <
typename TypeTag>
260 if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
265 this->linSys_.recoverSolutionWell(x, xw);
266 updateWellState(simulator, xw, well_state, deferred_logger);
273 template <
typename TypeTag>
278 std::vector<Scalar>& well_potentials,
281 const auto [compute_potential, bhp_controlled_well] =
284 if (!compute_potential) {
288 debug_cost_counter_ = 0;
289 bool converged_implicit =
false;
290 if (this->param_.local_well_solver_control_switching_) {
291 converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger);
292 if (!converged_implicit) {
293 deferred_logger.debug(
"Implicit potential calculations failed for well "
294 + this->name() +
", reverting to original aproach.");
297 if (!converged_implicit) {
299 const auto& summaryState = simulator.vanguard().summaryState();
300 if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
301 computeWellRatesAtBhpLimit(simulator, well_potentials, deferred_logger);
303 well_potentials = computeWellPotentialWithTHP(
304 well_state, simulator, deferred_logger);
307 deferred_logger.debug(
"Cost in iterations of finding well potential for well "
308 + this->name() +
": " + std::to_string(debug_cost_counter_));
310 this->checkNegativeWellPotentials(well_potentials,
311 this->param_.check_well_operability_,
318 template<
typename TypeTag>
322 std::vector<Scalar>& well_flux,
325 if (this->well_ecl_.isInjector()) {
326 const auto controls = this->well_ecl_.injectionControls(simulator.
vanguard().summaryState());
327 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
329 const auto controls = this->well_ecl_.productionControls(simulator.
vanguard().summaryState());
330 computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger);
334 template<
typename TypeTag>
336 MultisegmentWell<TypeTag>::
337 computeWellRatesWithBhp(
const Simulator& simulator,
339 std::vector<Scalar>& well_flux,
340 DeferredLogger& deferred_logger)
const
342 const int np = this->number_of_phases_;
344 well_flux.resize(np, 0.0);
345 const bool allow_cf = this->getAllowCrossFlow();
346 const int nseg = this->numberOfSegments();
347 const WellState<Scalar>& well_state = simulator.problem().wellModel().wellState();
348 const auto& ws = well_state.well(this->indexOfWell());
349 auto segments_copy = ws.segments;
350 segments_copy.scale_pressure(bhp);
351 const auto& segment_pressure = segments_copy.pressure;
352 for (
int seg = 0; seg < nseg; ++seg) {
353 for (
const int perf : this->segments_.perforations()[seg]) {
354 const int local_perf_index = this->pw_info_.activeToLocal(perf);
355 if (local_perf_index < 0)
357 const int cell_idx = this->well_cells_[local_perf_index];
358 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
360 std::vector<Scalar> mob(this->num_components_, 0.);
361 getMobility(simulator, local_perf_index, mob, deferred_logger);
362 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
363 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
364 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, intQuants, trans_mult, wellstate_nupcol);
365 const Scalar seg_pressure = segment_pressure[seg];
366 std::vector<Scalar> cq_s(this->num_components_, 0.);
367 Scalar perf_press = 0.0;
368 PerforationRates<Scalar> perf_rates;
369 computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
370 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
372 for(
int p = 0; p < np; ++p) {
373 well_flux[this->modelCompIdxToFlowCompIdx(p)] += cq_s[p];
377 this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
381 template<
typename TypeTag>
383 MultisegmentWell<TypeTag>::
384 computeWellRatesWithBhpIterations(
const Simulator& simulator,
386 std::vector<Scalar>& well_flux,
387 DeferredLogger& deferred_logger)
const
392 MultisegmentWell<TypeTag> well_copy(*
this);
393 well_copy.resetDampening();
395 well_copy.debug_cost_counter_ = 0;
398 WellState<Scalar> well_state_copy = simulator.problem().wellModel().wellState();
399 const auto& group_state = simulator.problem().wellModel().groupState();
400 auto& ws = well_state_copy.well(this->index_of_well_);
403 const auto& summary_state = simulator.vanguard().summaryState();
404 auto inj_controls = well_copy.well_ecl_.isInjector()
405 ? well_copy.well_ecl_.injectionControls(summary_state)
406 : Well::InjectionControls(0);
407 auto prod_controls = well_copy.well_ecl_.isProducer()
408 ? well_copy.well_ecl_.productionControls(summary_state) :
409 Well::ProductionControls(0);
412 if (well_copy.well_ecl_.isInjector()) {
413 inj_controls.bhp_limit = bhp;
414 ws.injection_cmode = Well::InjectorCMode::BHP;
416 prod_controls.bhp_limit = bhp;
417 ws.production_cmode = Well::ProducerCMode::BHP;
420 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
423 const int np = this->number_of_phases_;
425 for (
int phase = 0; phase < np; ++phase){
426 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
429 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
430 for (
int phase = 0; phase < np; ++phase) {
431 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
434 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
435 this->segments_.perforations(),
438 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
439 const double dt = simulator.timeStepSize();
441 well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
446 well_flux.resize(np, 0.0);
447 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
448 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
449 well_flux[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
451 debug_cost_counter_ += well_copy.debug_cost_counter_;
456 template<
typename TypeTag>
457 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
458 MultisegmentWell<TypeTag>::
459 computeWellPotentialWithTHP(
const WellState<Scalar>& well_state,
460 const Simulator& simulator,
461 DeferredLogger& deferred_logger)
const
463 std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
464 const auto& summary_state = simulator.vanguard().summaryState();
466 const auto& well = this->well_ecl_;
467 if (well.isInjector()) {
468 auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
469 if (bhp_at_thp_limit) {
470 const auto& controls = well.injectionControls(summary_state);
471 const Scalar bhp = std::min(*bhp_at_thp_limit,
472 static_cast<Scalar
>(controls.bhp_limit));
473 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
474 deferred_logger.debug(
"Converged thp based potential calculation for well "
475 + this->name() +
", at bhp = " + std::to_string(bhp));
477 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
478 "Failed in getting converged thp based potential calculation for well "
479 + this->name() +
". Instead the bhp based value is used");
480 const auto& controls = well.injectionControls(summary_state);
481 const Scalar bhp = controls.bhp_limit;
482 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
485 auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
486 well_state, simulator, summary_state, deferred_logger);
487 if (bhp_at_thp_limit) {
488 const auto& controls = well.productionControls(summary_state);
489 const Scalar bhp = std::max(*bhp_at_thp_limit,
490 static_cast<Scalar
>(controls.bhp_limit));
491 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
492 deferred_logger.debug(
"Converged thp based potential calculation for well "
493 + this->name() +
", at bhp = " + std::to_string(bhp));
495 deferred_logger.warning(
"FAILURE_GETTING_CONVERGED_POTENTIAL",
496 "Failed in getting converged thp based potential calculation for well "
497 + this->name() +
". Instead the bhp based value is used");
498 const auto& controls = well.productionControls(summary_state);
499 const Scalar bhp = controls.bhp_limit;
500 computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
507 template<
typename TypeTag>
509 MultisegmentWell<TypeTag>::
510 computeWellPotentialsImplicit(
const Simulator& simulator,
511 const WellState<Scalar>& well_state,
512 std::vector<Scalar>& well_potentials,
513 DeferredLogger& deferred_logger)
const
518 MultisegmentWell<TypeTag> well_copy(*
this);
519 well_copy.debug_cost_counter_ = 0;
522 WellState<Scalar> well_state_copy = well_state;
523 const auto& group_state = simulator.problem().wellModel().groupState();
524 auto& ws = well_state_copy.well(this->index_of_well_);
527 const auto& summary_state = simulator.vanguard().summaryState();
528 auto inj_controls = well_copy.well_ecl_.isInjector()
529 ? well_copy.well_ecl_.injectionControls(summary_state)
530 : Well::InjectionControls(0);
531 auto prod_controls = well_copy.well_ecl_.isProducer()
532 ? well_copy.well_ecl_.productionControls(summary_state)
533 : Well::ProductionControls(0);
536 well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls);
538 well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
541 const int np = this->number_of_phases_;
543 for (
int phase = 0; phase < np; ++phase){
544 trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
547 const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
548 for (
int phase = 0; phase < np; ++phase) {
549 ws.surface_rates[phase] = sign * ws.well_potentials[phase];
552 well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
553 this->segments_.perforations(),
556 well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
557 const double dt = simulator.timeStepSize();
559 bool converged =
false;
560 if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
561 converged = well_copy.solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
563 converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
567 well_potentials.clear();
568 well_potentials.resize(np, 0.0);
569 for (
int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
570 const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
571 well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
573 debug_cost_counter_ += well_copy.debug_cost_counter_;
577 template <
typename TypeTag>
579 MultisegmentWell<TypeTag>::
580 solveEqAndUpdateWellState(
const Simulator& simulator,
581 WellState<Scalar>& well_state,
582 DeferredLogger& deferred_logger)
584 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
589 const BVectorWell dx_well = this->linSys_.solve();
591 updateWellState(simulator, dx_well, well_state, deferred_logger);
593 catch(
const NumericalProblem& exp) {
597 deferred_logger.problem(
"In MultisegmentWell::solveEqAndUpdateWellState for well "
598 + this->name() +
": "+exp.what());
607 template <
typename TypeTag>
609 MultisegmentWell<TypeTag>::
610 computePerfCellPressDiffs(
const Simulator& simulator)
614 for (
int local_perf_index = 0; local_perf_index < this->number_of_local_perforations_; ++local_perf_index) {
617 std::vector<Scalar> kr(this->number_of_phases_, 0.0);
618 std::vector<Scalar> density(this->number_of_phases_, 0.0);
620 const int cell_idx = this->well_cells_[local_perf_index];
621 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
622 const auto& fs = intQuants.fluidState();
627 if (pu.phase_used[Water]) {
628 const int water_pos = pu.phase_pos[Water];
629 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
630 sum_kr += kr[water_pos];
631 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
634 if (pu.phase_used[Oil]) {
635 const int oil_pos = pu.phase_pos[Oil];
636 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
637 sum_kr += kr[oil_pos];
638 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
641 if (pu.phase_used[Gas]) {
642 const int gas_pos = pu.phase_pos[Gas];
643 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
644 sum_kr += kr[gas_pos];
645 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
648 assert(sum_kr != 0.);
651 Scalar average_density = 0.;
652 for (
int p = 0; p < this->number_of_phases_; ++p) {
653 average_density += kr[p] * density[p];
655 average_density /= sum_kr;
657 this->cell_perforation_pressure_diffs_[local_perf_index] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[local_perf_index];
665 template <
typename TypeTag>
667 MultisegmentWell<TypeTag>::
668 computeInitialSegmentFluids(
const Simulator& simulator)
670 for (
int seg = 0; seg < this->numberOfSegments(); ++seg) {
672 const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
673 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
674 segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
683 template <
typename TypeTag>
685 MultisegmentWell<TypeTag>::
686 updateWellState(
const Simulator& simulator,
687 const BVectorWell& dwells,
688 WellState<Scalar>& well_state,
689 DeferredLogger& deferred_logger,
690 const Scalar relaxation_factor)
692 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
694 const Scalar dFLimit = this->param_.dwell_fraction_max_;
695 const Scalar max_pressure_change = this->param_.max_pressure_change_ms_wells_;
696 const bool stop_or_zero_rate_target =
697 this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
698 this->primary_variables_.updateNewton(dwells,
701 stop_or_zero_rate_target,
702 max_pressure_change);
704 const auto& summary_state = simulator.vanguard().summaryState();
705 this->primary_variables_.copyToWellState(*
this, getRefDensity(),
711 auto& ws = well_state.well(this->index_of_well_);
712 this->segments_.copyPhaseDensities(ws.pu, ws.segments);
715 Base::calculateReservoirRates(simulator.vanguard().eclState().runspec().co2Storage(), well_state.well(this->index_of_well_));
722 template <
typename TypeTag>
724 MultisegmentWell<TypeTag>::
725 calculateExplicitQuantities(
const Simulator& simulator,
726 const WellState<Scalar>& well_state,
727 DeferredLogger& deferred_logger)
729 updatePrimaryVariables(simulator, well_state, deferred_logger);
730 computePerfCellPressDiffs(simulator);
731 computeInitialSegmentFluids(simulator);
738 template<
typename TypeTag>
740 MultisegmentWell<TypeTag>::
741 updateProductivityIndex(
const Simulator& simulator,
742 const WellProdIndexCalculator<Scalar>& wellPICalc,
743 WellState<Scalar>& well_state,
744 DeferredLogger& deferred_logger)
const
746 auto fluidState = [&simulator,
this](
const int local_perf_index)
748 const auto cell_idx = this->well_cells_[local_perf_index];
749 return simulator.model()
750 .intensiveQuantities(cell_idx, 0).fluidState();
753 const int np = this->number_of_phases_;
754 auto setToZero = [np](Scalar* x) ->
void
756 std::fill_n(x, np, 0.0);
759 auto addVector = [np](
const Scalar* src, Scalar* dest) ->
void
761 std::transform(src, src + np, dest, dest, std::plus<>{});
764 auto& ws = well_state.well(this->index_of_well_);
765 auto& perf_data = ws.perf_data;
766 auto* connPI = perf_data.prod_index.data();
767 auto* wellPI = ws.productivity_index.data();
771 const auto preferred_phase = this->well_ecl_.getPreferredPhase();
772 auto subsetPerfID = 0;
774 for (
const auto& perf : *this->perf_data_){
775 auto allPerfID = perf.ecl_index;
777 auto connPICalc = [&wellPICalc, allPerfID](
const Scalar mobility) -> Scalar
779 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
782 std::vector<Scalar> mob(this->num_components_, 0.0);
787 getMobility(simulator,
static_cast<int>(subsetPerfID), mob, deferred_logger);
789 const auto& fs = fluidState(subsetPerfID);
792 if (this->isInjector()) {
793 this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
794 mob, connPI, deferred_logger);
797 this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
800 addVector(connPI, wellPI);
807 const auto& comm = this->parallel_well_info_.communication();
808 if (comm.size() > 1) {
809 comm.sum(wellPI, np);
812 assert (
static_cast<int>(subsetPerfID) == this->number_of_local_perforations_ &&
813 "Internal logic error in processing connections for PI/II");
820 template<
typename TypeTag>
821 typename MultisegmentWell<TypeTag>::Scalar
822 MultisegmentWell<TypeTag>::
823 connectionDensity(
const int globalConnIdx,
824 [[maybe_unused]]
const int openConnIdx)
const
829 const auto segNum = this->wellEcl()
830 .getConnections()[globalConnIdx].segment();
832 const auto segIdx = this->wellEcl()
833 .getSegments().segmentNumberToIndex(segNum);
835 return this->segments_.density(segIdx).value();
842 template<
typename TypeTag>
844 MultisegmentWell<TypeTag>::
845 addWellContributions(SparseMatrixAdapter& jacobian)
const
847 if (this->number_of_local_perforations_ == 0) {
851 this->linSys_.extract(jacobian);
855 template<
typename TypeTag>
857 MultisegmentWell<TypeTag>::
858 addWellPressureEquations(PressureMatrix& jacobian,
859 const BVector& weights,
860 const int pressureVarIndex,
861 const bool use_well_weights,
862 const WellState<Scalar>& well_state)
const
864 if (this->number_of_local_perforations_ == 0) {
869 this->linSys_.extractCPRPressureMatrix(jacobian,
879 template<
typename TypeTag>
880 template<
class Value>
882 MultisegmentWell<TypeTag>::
883 computePerfRate(
const Value& pressure_cell,
886 const std::vector<Value>& b_perfcells,
887 const std::vector<Value>& mob_perfcells,
888 const std::vector<Scalar>& Tw,
890 const Value& segment_pressure,
891 const Value& segment_density,
892 const bool& allow_cf,
893 const std::vector<Value>& cmix_s,
894 std::vector<Value>& cq_s,
896 PerforationRates<Scalar>& perf_rates,
897 DeferredLogger& deferred_logger)
const
899 const int local_perf_index = this->pw_info_.activeToLocal(perf);
900 if (local_perf_index < 0)
904 const Value perf_seg_press_diff = this->gravity() * segment_density *
905 this->segments_.local_perforation_depth_diff(local_perf_index);
907 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
911 perf_press = segment_pressure + perf_seg_press_diff;
914 const Value cell_press_at_perf = pressure_cell - cell_perf_press_diff;
917 const Value drawdown = cell_press_at_perf - perf_press;
920 if (drawdown > 0.0) {
922 if (!allow_cf && this->isInjector()) {
927 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
928 const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
929 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
932 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
933 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
934 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
935 const Value cq_s_oil = cq_s[oilCompIdx];
936 const Value cq_s_gas = cq_s[gasCompIdx];
937 cq_s[gasCompIdx] += rs * cq_s_oil;
938 cq_s[oilCompIdx] += rv * cq_s_gas;
942 if (!allow_cf && this->isProducer()) {
947 Value total_mob = mob_perfcells[0];
948 for (
int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
949 total_mob += mob_perfcells[comp_idx];
953 Value volume_ratio = 0.0;
954 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
955 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
956 volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
959 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
960 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
961 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
966 const Value d = 1.0 - rv * rs;
968 if (getValue(d) == 0.0) {
969 OPM_DEFLOG_PROBLEM(NumericalProblem,
970 fmt::format(
"Zero d value obtained for well {} "
971 "during flux calculation with rs {} and rv {}",
972 this->name(), rs, rv),
976 const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
977 volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
979 const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
980 volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
982 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
983 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
984 volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
986 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
987 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
988 volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
992 for (
int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
993 const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown);
994 Value cqt_is = cqt_i / volume_ratio;
995 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
1000 if (this->isProducer()) {
1001 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1002 const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1003 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1012 const Scalar d = 1.0 - getValue(rv) * getValue(rs);
1015 perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
1018 perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
1023 template <
typename TypeTag>
1024 template<
class Value>
1026 MultisegmentWell<TypeTag>::
1027 computePerfRate(
const IntensiveQuantities& int_quants,
1028 const std::vector<Value>& mob_perfcells,
1029 const std::vector<Scalar>& Tw,
1032 const Value& segment_pressure,
1033 const bool& allow_cf,
1034 std::vector<Value>& cq_s,
1036 PerforationRates<Scalar>& perf_rates,
1037 DeferredLogger& deferred_logger)
const
1040 auto obtain = [
this](
const Eval& value)
1042 if constexpr (std::is_same_v<Value, Scalar>) {
1043 static_cast<void>(
this);
1044 return getValue(value);
1046 return this->extendEval(value);
1049 auto obtainN = [](
const auto& value)
1051 if constexpr (std::is_same_v<Value, Scalar>) {
1052 return getValue(value);
1057 const auto& fs = int_quants.fluidState();
1059 const Value pressure_cell = obtain(this->getPerfCellPressure(fs));
1060 const Value rs = obtain(fs.Rs());
1061 const Value rv = obtain(fs.Rv());
1064 std::vector<Value> b_perfcells(this->num_components_, 0.0);
1066 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1067 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1071 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1072 b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
1075 std::vector<Value> cmix_s(this->numComponents(), 0.0);
1076 for (
int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
1077 cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
1080 this->computePerfRate(pressure_cell,
1088 obtainN(this->segments_.density(seg)),
1097 template <
typename TypeTag>
1099 MultisegmentWell<TypeTag>::
1100 computeSegmentFluidProperties(
const Simulator& simulator, DeferredLogger& deferred_logger)
1109 EvalWell temperature;
1110 EvalWell saltConcentration;
1116 int pvt_region_index;
1118 Scalar fsTemperature;
1119 using SaltConcType =
typename std::decay<
decltype(std::declval<
decltype(simulator.model().intensiveQuantities(0, 0).fluidState())>().saltConcentration())>::type;
1120 SaltConcType fsSaltConcentration;
1122 if (this->well_cells_.size() > 0) {
1126 const int cell_idx = this->well_cells_[0];
1127 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1128 const auto& fs = intQuants.fluidState();
1130 fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1131 fsSaltConcentration = fs.saltConcentration();
1132 pvt_region_index = fs.pvtRegionIndex();
1137 fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
1138 temperature.setValue(fsTemperature);
1140 fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
1141 saltConcentration = this->extendEval(fsSaltConcentration);
1143 pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
1145 this->segments_.computeFluidProperties(temperature,
1147 this->primary_variables_,
1152 template <
typename TypeTag>
1153 template<
class Value>
1155 MultisegmentWell<TypeTag>::
1156 getMobility(
const Simulator& simulator,
1157 const int local_perf_index,
1158 std::vector<Value>& mob,
1159 DeferredLogger& deferred_logger)
const
1161 auto obtain = [
this](
const Eval& value)
1163 if constexpr (std::is_same_v<Value, Scalar>) {
1164 static_cast<void>(
this);
1165 return getValue(value);
1167 return this->extendEval(value);
1171 WellInterface<TypeTag>::getMobility(simulator, local_perf_index, mob, obtain, deferred_logger);
1173 if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
1174 const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index;
1175 const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
1176 const int seg = this->segmentNumberToIndex(con.segment());
1180 const Scalar segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
1181 const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
1182 * this->segments_.local_perforation_depth_diff(local_perf_index);
1183 const Scalar perf_press = segment_pres + perf_seg_press_diff;
1184 const Scalar multiplier = this->getInjMult(local_perf_index, segment_pres, perf_press, deferred_logger);
1185 for (std::size_t i = 0; i < mob.size(); ++i) {
1186 mob[i] *= multiplier;
1193 template<
typename TypeTag>
1194 typename MultisegmentWell<TypeTag>::Scalar
1195 MultisegmentWell<TypeTag>::
1196 getRefDensity()
const
1198 return this->segments_.getRefDensity();
1201 template<
typename TypeTag>
1203 MultisegmentWell<TypeTag>::
1204 checkOperabilityUnderBHPLimit(
const WellState<Scalar>& ,
1205 const Simulator& simulator,
1206 DeferredLogger& deferred_logger)
1208 const auto& summaryState = simulator.vanguard().summaryState();
1209 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1212 const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1213 if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1216 Scalar total_ipr_mass_rate = 0.0;
1217 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1219 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1223 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1224 const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1226 const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1227 total_ipr_mass_rate += ipr_rate * rho;
1229 if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1230 this->operability_status_.operable_under_only_bhp_limit =
false;
1234 if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1238 std::vector<Scalar> well_rates_bhp_limit;
1239 computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1241 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1242 const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
1244 this->getRefDensity(),
1245 this->wellEcl().alq_value(summaryState),
1248 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1249 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1260 this->operability_status_.operable_under_only_bhp_limit =
true;
1261 this->operability_status_.obey_thp_limit_under_bhp_limit =
false;
1267 template<
typename TypeTag>
1269 MultisegmentWell<TypeTag>::
1270 updateIPR(
const Simulator& simulator, DeferredLogger& deferred_logger)
const
1275 std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1276 std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1278 const int nseg = this->numberOfSegments();
1279 std::vector<Scalar> seg_dp(nseg, 0.0);
1280 for (
int seg = 0; seg < nseg; ++seg) {
1282 const Scalar dp = this->getSegmentDp(seg,
1283 this->segments_.density(seg).value(),
1286 for (
const int perf : this->segments_.perforations()[seg]) {
1287 const int local_perf_index = this->pw_info_.activeToLocal(perf);
1288 if (local_perf_index < 0)
1290 std::vector<Scalar> mob(this->num_components_, 0.0);
1293 getMobility(simulator, local_perf_index, mob, deferred_logger);
1295 const int cell_idx = this->well_cells_[local_perf_index];
1296 const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, 0);
1297 const auto& fs = int_quantities.fluidState();
1299 const Scalar perf_seg_press_diff = this->segments_.getPressureDiffSegLocalPerf(seg, local_perf_index);
1301 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
1302 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1305 std::vector<Scalar> b_perf(this->num_components_);
1306 for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1307 if (!FluidSystem::phaseIsActive(phase)) {
1310 const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1311 b_perf[comp_idx] = fs.invB(phase).value();
1315 const Scalar h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
1316 const Scalar pressure_diff = pressure_cell - h_perf;
1319 if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1320 deferred_logger.debug(
"CROSSFLOW_IPR",
1321 "cross flow found when updateIPR for well " + this->name());
1325 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
1326 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1327 const std::vector<Scalar> tw_perf = this->wellIndex(local_perf_index, int_quantities, trans_mult, wellstate_nupcol);
1328 std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
1329 std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
1330 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1331 const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
1332 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1333 ipr_b_perf[comp_idx] += tw_mob;
1337 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1338 const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1339 const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1340 const Scalar rs = (fs.Rs()).value();
1341 const Scalar rv = (fs.Rv()).value();
1343 const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1344 const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1346 ipr_a_perf[gas_comp_idx] += dis_gas_a;
1347 ipr_a_perf[oil_comp_idx] += vap_oil_a;
1349 const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1350 const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1352 ipr_b_perf[gas_comp_idx] += dis_gas_b;
1353 ipr_b_perf[oil_comp_idx] += vap_oil_b;
1356 for (std::size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1357 this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1358 this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1362 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1363 this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1366 template<
typename TypeTag>
1368 MultisegmentWell<TypeTag>::
1369 updateIPRImplicit(
const Simulator& simulator,
1370 WellState<Scalar>& well_state,
1371 DeferredLogger& deferred_logger)
1380 auto rates = well_state.well(this->index_of_well_).surface_rates;
1382 for (std::size_t p = 0; p < rates.size(); ++p) {
1383 zero_rates &= rates[p] == 0.0;
1385 auto& ws = well_state.well(this->index_of_well_);
1387 const auto msg = fmt::format(
"updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name());
1388 deferred_logger.debug(msg);
1400 const auto& group_state = simulator.problem().wellModel().groupState();
1402 std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
1403 std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
1405 auto inj_controls = Well::InjectionControls(0);
1406 auto prod_controls = Well::ProductionControls(0);
1407 prod_controls.addControl(Well::ProducerCMode::BHP);
1408 prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
1411 const auto cmode = ws.production_cmode;
1412 ws.production_cmode = Well::ProducerCMode::BHP;
1413 const double dt = simulator.timeStepSize();
1414 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1416 BVectorWell rhs(this->numberOfSegments());
1418 rhs[0][SPres] = -1.0;
1420 const BVectorWell x_well = this->linSys_.solve(rhs);
1421 constexpr int num_eq = MSWEval::numWellEq;
1422 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
1423 const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
1424 const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
1425 for (
size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
1427 ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
1429 ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
1432 ws.production_cmode = cmode;
1435 template<
typename TypeTag>
1437 MultisegmentWell<TypeTag>::
1438 checkOperabilityUnderTHPLimit(
const Simulator& simulator,
1439 const WellState<Scalar>& well_state,
1440 DeferredLogger& deferred_logger)
1442 const auto& summaryState = simulator.vanguard().summaryState();
1443 const auto obtain_bhp = this->isProducer()
1444 ? computeBhpAtThpLimitProd(
1445 well_state, simulator, summaryState, deferred_logger)
1446 : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger);
1449 this->operability_status_.can_obtain_bhp_with_thp_limit =
true;
1451 const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
1452 this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1454 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1455 if (this->isProducer() && *obtain_bhp < thp_limit) {
1456 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1457 +
" bars is SMALLER than thp limit "
1458 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1459 +
" bars as a producer for well " + this->name();
1460 deferred_logger.debug(msg);
1462 else if (this->isInjector() && *obtain_bhp > thp_limit) {
1463 const std::string msg =
" obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1464 +
" bars is LARGER than thp limit "
1465 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1466 +
" bars as a injector for well " + this->name();
1467 deferred_logger.debug(msg);
1472 this->operability_status_.can_obtain_bhp_with_thp_limit =
false;
1473 this->operability_status_.obey_bhp_limit_with_thp_limit =
false;
1474 if (!this->wellIsStopped()) {
1475 const Scalar thp_limit = this->getTHPConstraint(summaryState);
1476 deferred_logger.debug(
" could not find bhp value at thp limit "
1477 + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1478 +
" bar for well " + this->name() +
", the well might need to be closed ");
1487 template<
typename TypeTag>
1489 MultisegmentWell<TypeTag>::
1490 iterateWellEqWithControl(
const Simulator& simulator,
1492 const Well::InjectionControls& inj_controls,
1493 const Well::ProductionControls& prod_controls,
1494 WellState<Scalar>& well_state,
1495 const GroupState<Scalar>& group_state,
1496 DeferredLogger& deferred_logger)
1498 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return true;
1500 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1504 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1509 std::vector<std::vector<Scalar> > residual_history;
1510 std::vector<Scalar> measure_history;
1513 Scalar relaxation_factor = 1.;
1514 bool converged =
false;
1515 bool relax_convergence =
false;
1516 this->regularize_ =
false;
1517 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1519 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1521 BVectorWell dx_well;
1523 dx_well = this->linSys_.solve();
1525 catch(
const NumericalProblem& exp) {
1529 deferred_logger.problem(
"In MultisegmentWell::iterateWellEqWithControl for well "
1530 + this->name() +
": "+exp.what());
1534 if (it > this->param_.strict_inner_iter_wells_) {
1535 relax_convergence =
true;
1536 this->regularize_ =
true;
1539 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1540 if (report.converged()) {
1547 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1551 residual_history.push_back(residuals);
1552 measure_history.push_back(this->getResidualMeasureValue(well_state,
1553 residual_history[it],
1554 this->param_.tolerance_wells_,
1555 this->param_.tolerance_pressure_ms_wells_,
1558 bool min_relaxation_reached = this->update_relaxation_factor(measure_history, relaxation_factor, this->regularize_, deferred_logger);
1559 if (min_relaxation_reached || this->repeatedStagnation(measure_history, this->regularize_, deferred_logger)) {
1561 const auto reportStag = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger,
true);
1562 if (reportStag.converged()) {
1564 std::string message = fmt::format(
"Well stagnates/oscillates but {} manages to get converged with relaxed tolerances in {} inner iterations."
1566 deferred_logger.debug(message);
1572 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1577 std::ostringstream sstr;
1578 sstr <<
" Well " << this->name() <<
" converged in " << it <<
" inner iterations.";
1579 if (relax_convergence)
1580 sstr <<
" (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ <<
" iterations)";
1581 deferred_logger.debug(sstr.str());
1583 std::ostringstream sstr;
1584 sstr <<
" Well " << this->name() <<
" did not converge in " << it <<
" inner iterations.";
1585#define EXTRA_DEBUG_MSW 0
1587 sstr <<
"***** Outputting the residual history for well " << this->name() <<
" during inner iterations:";
1588 for (
int i = 0; i < it; ++i) {
1589 const auto& residual = residual_history[i];
1590 sstr <<
" residual at " << i <<
"th iteration ";
1591 for (
const auto& res : residual) {
1594 sstr <<
" " << measure_history[i] <<
" \n";
1597#undef EXTRA_DEBUG_MSW
1598 deferred_logger.debug(sstr.str());
1605 template<
typename TypeTag>
1607 MultisegmentWell<TypeTag>::
1608 iterateWellEqWithSwitching(
const Simulator& simulator,
1610 const Well::InjectionControls& inj_controls,
1611 const Well::ProductionControls& prod_controls,
1612 WellState<Scalar>& well_state,
1613 const GroupState<Scalar>& group_state,
1614 DeferredLogger& deferred_logger,
1615 const bool fixed_control ,
1616 const bool fixed_status )
1618 const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1622 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1627 std::vector<std::vector<Scalar> > residual_history;
1628 std::vector<Scalar> measure_history;
1631 Scalar relaxation_factor = 1.;
1632 bool converged =
false;
1633 bool relax_convergence =
false;
1634 this->regularize_ =
false;
1635 const auto& summary_state = simulator.vanguard().summaryState();
1640 const int min_its_after_switch = 3;
1641 int its_since_last_switch = min_its_after_switch;
1642 int switch_count= 0;
1644 const auto well_status_orig = this->wellStatus_;
1645 const auto operability_orig = this->operability_status_;
1647 const bool allow_open = this->well_ecl_.getStatus() == WellStatus::OPEN &&
1648 well_state.well(this->index_of_well_).status == WellStatus::OPEN;
1650 const bool allow_switching = !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
1651 (!fixed_control || !fixed_status) && allow_open;
1652 bool final_check =
false;
1654 this->operability_status_.resetOperability();
1655 this->operability_status_.solvable =
true;
1657 for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1658 ++its_since_last_switch;
1659 if (allow_switching && its_since_last_switch >= min_its_after_switch){
1660 const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
1661 bool changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
1662 inj_controls, prod_controls, wqTotal,
1663 deferred_logger, fixed_control,
1666 its_since_last_switch = 0;
1669 if (!changed && final_check) {
1672 final_check =
false;
1676 assembleWellEqWithoutIteration(simulator, dt, inj_controls, prod_controls,
1677 well_state, group_state, deferred_logger);
1679 const BVectorWell dx_well = this->linSys_.solve();
1681 if (it > this->param_.strict_inner_iter_wells_) {
1682 relax_convergence =
true;
1683 this->regularize_ =
true;
1686 const auto report = getWellConvergence(simulator, well_state, Base::B_avg_, deferred_logger, relax_convergence);
1687 converged = report.converged();
1688 if (this->parallel_well_info_.communication().size() > 1 && converged != this->parallel_well_info_.communication().min(converged)) {
1689 OPM_THROW(std::runtime_error, fmt::format(
"Misalignment of the parallel simulation run in iterateWellEqWithSwitching - the well calculation succeeded on rank {} but failed on other ranks.", this->parallel_well_info_.communication().rank()));
1694 if (switch_count > 0 && its_since_last_switch < min_its_after_switch) {
1696 its_since_last_switch = min_its_after_switch;
1704 const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
1708 residual_history.push_back(residuals);
1712 measure_history.push_back(this->getResidualMeasureValue(well_state,
1713 residual_history[it],
1714 this->param_.tolerance_wells_,
1715 this->param_.tolerance_pressure_ms_wells_,
1717 bool min_relaxation_reached = this->update_relaxation_factor(measure_history, relaxation_factor, this->regularize_, deferred_logger);
1718 if (min_relaxation_reached || this->repeatedStagnation(measure_history, this->regularize_, deferred_logger)) {
1723 updateWellState(simulator, dx_well, well_state, deferred_logger, relaxation_factor);
1727 if (allow_switching){
1729 const bool is_stopped = this->wellIsStopped();
1730 if (this->wellHasTHPConstraints(summary_state)){
1731 this->operability_status_.can_obtain_bhp_with_thp_limit = !is_stopped;
1732 this->operability_status_.obey_thp_limit_under_bhp_limit = !is_stopped;
1734 this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
1737 std::string message = fmt::format(
" Well {} converged in {} inner iterations ("
1738 "{} control/status switches).", this->name(), it, switch_count);
1739 if (relax_convergence) {
1740 message.append(fmt::format(
" (A relaxed tolerance was used after {} iterations)",
1741 this->param_.strict_inner_iter_wells_));
1743 deferred_logger.debug(message);
1745 this->wellStatus_ = well_status_orig;
1746 this->operability_status_ = operability_orig;
1747 const std::string message = fmt::format(
" Well {} did not converge in {} inner iterations ("
1748 "{} control/status switches).", this->name(), it, switch_count);
1749 deferred_logger.debug(message);
1750 this->primary_variables_.outputLowLimitPressureSegments(deferred_logger);
1757 template<
typename TypeTag>
1759 MultisegmentWell<TypeTag>::
1760 assembleWellEqWithoutIteration(
const Simulator& simulator,
1762 const Well::InjectionControls& inj_controls,
1763 const Well::ProductionControls& prod_controls,
1764 WellState<Scalar>& well_state,
1765 const GroupState<Scalar>& group_state,
1766 DeferredLogger& deferred_logger)
1768 if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
1771 this->segments_.updateUpwindingSegments(this->primary_variables_);
1774 computeSegmentFluidProperties(simulator, deferred_logger);
1777 this->linSys_.clear();
1779 auto& ws = well_state.well(this->index_of_well_);
1780 ws.phase_mixing_rates.fill(0.0);
1787 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
1789 const int nseg = this->numberOfSegments();
1791 for (
int seg = 0; seg < nseg; ++seg) {
1793 const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
1794 auto& perf_data = ws.perf_data;
1795 auto& perf_rates = perf_data.phase_rates;
1796 auto& perf_press_state = perf_data.pressure;
1797 for (
const int perf : this->segments_.perforations()[seg]) {
1798 const int local_perf_index = this->pw_info_.activeToLocal(perf);
1799 if (local_perf_index < 0)
1801 const int cell_idx = this->well_cells_[local_perf_index];
1802 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
1803 std::vector<EvalWell> mob(this->num_components_, 0.0);
1804 getMobility(simulator, local_perf_index, mob, deferred_logger);
1805 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
1806 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
1807 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
1808 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1809 EvalWell perf_press;
1810 PerforationRates<Scalar> perfRates;
1811 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
1812 allow_cf, cq_s, perf_press, perfRates, deferred_logger);
1815 if (this->isProducer()) {
1816 ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas;
1817 ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil;
1818 perf_data.phase_mixing_rates[local_perf_index][ws.dissolved_gas] = perfRates.dis_gas;
1819 perf_data.phase_mixing_rates[local_perf_index][ws.vaporized_oil] = perfRates.vap_oil;
1823 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1824 perf_rates[local_perf_index*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1826 perf_press_state[local_perf_index] = perf_press.value();
1828 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1830 const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1832 this->connectionRates_[local_perf_index][comp_idx] = Base::restrictEval(cq_s_effective);
1834 MultisegmentWellAssemble(*this).
1835 assemblePerforationEq(seg, local_perf_index, comp_idx, cq_s_effective, this->linSys_);
1841 const auto& comm = this->parallel_well_info_.communication();
1842 comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
1845 if (this->parallel_well_info_.communication().size() > 1) {
1847 this->linSys_.sumDistributed(this->parallel_well_info_.communication());
1849 for (
int seg = 0; seg < nseg; ++seg) {
1853 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg);
1858 const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
1860 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1861 const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
1862 - segment_fluid_initial_[seg][comp_idx]) / dt;
1863 MultisegmentWellAssemble(*this).
1864 assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
1869 const int seg_upwind = this->segments_.upwinding_segment(seg);
1870 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1871 const EvalWell segment_rate =
1872 this->primary_variables_.getSegmentRateUpwinding(seg,
1875 this->well_efficiency_factor_;
1876 MultisegmentWellAssemble(*this).
1877 assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
1883 for (
const int inlet : this->segments_.inlets()[seg]) {
1884 const int inlet_upwind = this->segments_.upwinding_segment(inlet);
1885 for (
int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1886 const EvalWell inlet_rate =
1887 this->primary_variables_.getSegmentRateUpwinding(inlet,
1890 this->well_efficiency_factor_;
1891 MultisegmentWellAssemble(*this).
1892 assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
1899 const auto& summaryState = simulator.vanguard().summaryState();
1900 const Schedule& schedule = simulator.vanguard().schedule();
1901 const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
1902 MultisegmentWellAssemble(*this).
1903 assembleControlEq(well_state,
1910 this->primary_variables_,
1912 stopped_or_zero_target,
1915 const UnitSystem& unit_system = simulator.vanguard().eclState().getDeckUnitSystem();
1916 const auto& summary_state = simulator.vanguard().summaryState();
1917 this->assemblePressureEq(seg, unit_system, well_state, summary_state, this->param_.use_average_density_ms_wells_, deferred_logger);
1921 this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1922 this->linSys_.createSolver();
1928 template<
typename TypeTag>
1930 MultisegmentWell<TypeTag>::
1931 openCrossFlowAvoidSingularity(
const Simulator& simulator)
const
1933 return !this->getAllowCrossFlow() && allDrawDownWrongDirection(simulator);
1937 template<
typename TypeTag>
1939 MultisegmentWell<TypeTag>::
1940 allDrawDownWrongDirection(
const Simulator& simulator)
const
1942 bool all_drawdown_wrong_direction =
true;
1943 const int nseg = this->numberOfSegments();
1945 for (
int seg = 0; seg < nseg; ++seg) {
1946 const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
1947 for (
const int perf : this->segments_.perforations()[seg]) {
1948 const int local_perf_index = this->pw_info_.activeToLocal(perf);
1949 if (local_perf_index < 0)
1952 const int cell_idx = this->well_cells_[local_perf_index];
1953 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
1954 const auto& fs = intQuants.fluidState();
1957 const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegLocalPerf(seg, local_perf_index);
1959 const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index];
1961 const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
1962 const Scalar perf_press = pressure_cell - cell_perf_press_diff;
1965 const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1970 if ( (drawdown < 0. && this->isInjector()) ||
1971 (drawdown > 0. && this->isProducer()) ) {
1972 all_drawdown_wrong_direction =
false;
1977 const auto& comm = this->parallel_well_info_.communication();
1978 if (comm.size() > 1)
1980 all_drawdown_wrong_direction =
1981 (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1984 return all_drawdown_wrong_direction;
1990 template<
typename TypeTag>
1992 MultisegmentWell<TypeTag>::
1993 updateWaterThroughput(
const double , WellState<Scalar>& )
const
2001 template<
typename TypeTag>
2002 typename MultisegmentWell<TypeTag>::EvalWell
2003 MultisegmentWell<TypeTag>::
2004 getSegmentSurfaceVolume(
const Simulator& simulator,
const int seg_idx)
const
2006 EvalWell temperature;
2007 EvalWell saltConcentration;
2008 int pvt_region_index;
2010 Scalar fsTemperature;
2011 using SaltConcType =
typename std::decay<
decltype(std::declval<
decltype(simulator.model().intensiveQuantities(0, 0).fluidState())>().saltConcentration())>::type;
2012 SaltConcType fsSaltConcentration;
2014 if (this->well_cells_.size() > 0) {
2018 const int cell_idx = this->well_cells_[0];
2019 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, 0);
2020 const auto& fs = intQuants.fluidState();
2022 fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
2023 fsSaltConcentration = fs.saltConcentration();
2024 pvt_region_index = fs.pvtRegionIndex();
2029 fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
2030 temperature.setValue(fsTemperature);
2032 fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
2033 saltConcentration = this->extendEval(fsSaltConcentration);
2035 pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
2037 return this->segments_.getSurfaceVolume(temperature,
2039 this->primary_variables_,
2045 template<
typename TypeTag>
2046 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2047 MultisegmentWell<TypeTag>::
2048 computeBhpAtThpLimitProd(
const WellState<Scalar>& well_state,
2049 const Simulator& simulator,
2050 const SummaryState& summary_state,
2051 DeferredLogger& deferred_logger)
const
2053 return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
2056 this->getALQ(well_state),
2063 template<
typename TypeTag>
2064 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2065 MultisegmentWell<TypeTag>::
2066 computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
2067 const SummaryState& summary_state,
2068 const Scalar alq_value,
2069 DeferredLogger& deferred_logger,
2070 bool iterate_if_no_solution)
const
2074 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2080 std::vector<Scalar> rates(3);
2081 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2085 auto bhpAtLimit = WellBhpThpCalculator(*this).
2086 computeBhpAtThpLimitProd(frates,
2088 this->maxPerfPress(simulator),
2089 this->getRefDensity(),
2091 this->getTHPConstraint(summary_state),
2097 if (!iterate_if_no_solution)
2098 return std::nullopt;
2100 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2104 std::vector<Scalar> rates(3);
2105 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2109 return WellBhpThpCalculator(*this).
2110 computeBhpAtThpLimitProd(fratesIter,
2112 this->maxPerfPress(simulator),
2113 this->getRefDensity(),
2115 this->getTHPConstraint(summary_state),
2119 template<
typename TypeTag>
2120 std::optional<typename MultisegmentWell<TypeTag>::Scalar>
2121 MultisegmentWell<TypeTag>::
2122 computeBhpAtThpLimitInj(
const Simulator& simulator,
2123 const SummaryState& summary_state,
2124 DeferredLogger& deferred_logger)
const
2127 auto frates = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2133 std::vector<Scalar> rates(3);
2134 computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
2138 auto bhpAtLimit = WellBhpThpCalculator(*this).
2139 computeBhpAtThpLimitInj(frates,
2141 this->getRefDensity(),
2150 auto fratesIter = [
this, &simulator, &deferred_logger](
const Scalar bhp) {
2154 std::vector<Scalar> rates(3);
2155 computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
2159 return WellBhpThpCalculator(*this).
2160 computeBhpAtThpLimitInj(fratesIter,
2162 this->getRefDensity(),
2173 template<
typename TypeTag>
2174 typename MultisegmentWell<TypeTag>::Scalar
2175 MultisegmentWell<TypeTag>::
2176 maxPerfPress(
const Simulator& simulator)
const
2178 Scalar max_pressure = 0.0;
2179 const int nseg = this->numberOfSegments();
2180 for (
int seg = 0; seg < nseg; ++seg) {
2181 for (
const int perf : this->segments_.perforations()[seg]) {
2182 const int local_perf_index = this->pw_info_.activeToLocal(perf);
2183 if (local_perf_index < 0)
2186 const int cell_idx = this->well_cells_[local_perf_index];
2187 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2188 const auto& fs = int_quants.fluidState();
2189 Scalar pressure_cell = this->getPerfCellPressure(fs).value();
2190 max_pressure = std::max(max_pressure, pressure_cell);
2193 return max_pressure;
2200 template<
typename TypeTag>
2201 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2207 std::vector<Scalar> well_q_s(this->num_components_, 0.0);
2208 const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
2209 const int nseg = this->numberOfSegments();
2210 for (
int seg = 0; seg < nseg; ++seg) {
2212 const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
2213 for (
const int perf : this->segments_.perforations()[seg]) {
2214 const int local_perf_index = this->pw_info_.activeToLocal(perf);
2215 if (local_perf_index < 0)
2218 const int cell_idx = this->well_cells_[local_perf_index];
2219 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, 0);
2220 std::vector<Scalar> mob(this->num_components_, 0.0);
2221 getMobility(simulator, local_perf_index, mob, deferred_logger);
2222 const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
2223 const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
2224 const std::vector<Scalar> Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
2225 std::vector<Scalar> cq_s(this->num_components_, 0.0);
2226 Scalar perf_press = 0.0;
2228 computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
2229 allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
2230 for (
int comp = 0; comp < this->num_components_; ++comp) {
2231 well_q_s[comp] += cq_s[comp];
2235 const auto& comm = this->parallel_well_info_.communication();
2236 if (comm.size() > 1)
2238 comm.sum(well_q_s.data(), well_q_s.size());
2244 template <
typename TypeTag>
2245 std::vector<typename MultisegmentWell<TypeTag>::Scalar>
2249 const int num_seg = this->numberOfSegments();
2250 constexpr int num_eq = MSWEval::numWellEq;
2251 std::vector<Scalar> retval(num_seg * num_eq);
2252 for (
int ii = 0; ii < num_seg; ++ii) {
2253 const auto& pv = this->primary_variables_.value(ii);
2254 std::copy(pv.begin(), pv.end(), retval.begin() + ii * num_eq);
2262 template <
typename TypeTag>
2264 MultisegmentWell<TypeTag>::
2265 setPrimaryVars(
typename std::vector<Scalar>::const_iterator it)
2267 const int num_seg = this->numberOfSegments();
2268 constexpr int num_eq = MSWEval::numWellEq;
2269 std::array<Scalar, num_eq> tmp;
2270 for (
int ii = 0; ii < num_seg; ++ii) {
2271 const auto start = it + ii * num_eq;
2272 std::copy(start, start + num_eq, tmp.begin());
2273 this->primary_variables_.setValue(ii, tmp);
2275 return num_seg * num_eq;
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
Definition DeferredLogger.hpp:57
Definition GroupState.hpp:43
Definition MultisegmentWell.hpp:38
void apply(const BVector &x, BVector &Ax) const override
Ax = Ax - C D^-1 B x.
Definition MultisegmentWell_impl.hpp:220
Manages the initializing and running of time dependent problems.
Definition simulator.hh:97
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:281
Definition WellInterfaceGeneric.hpp:53
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:66
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition PerforationData.hpp:41