22 #include <opm/autodiff/MSWellHelpers.hpp> 28 template <
typename TypeTag>
29 MultisegmentWell<TypeTag>::
30 MultisegmentWell(
const Well* well,
const int time_step,
const Wells* wells,
const ModelParameters& param)
31 : Base(well, time_step, wells, param)
32 , segment_perforations_(numberOfSegments())
33 , segment_inlets_(numberOfSegments())
34 , cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
35 , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0)
36 , perforation_segment_depth_diffs_(number_of_perforations_, 0.0)
37 , segment_comp_initial_(numberOfSegments(),
std::vector<double>(numComponents(), 0.0))
38 , segment_densities_(numberOfSegments(), 0.0)
39 , segment_viscosities_(numberOfSegments(), 0.0)
40 , segment_mass_rates_(numberOfSegments(), 0.0)
41 , segment_depth_diffs_(numberOfSegments(), 0.0)
45 OPM_THROW(std::runtime_error,
"solvent is not supported by multisegment well yet");
49 OPM_THROW(std::runtime_error,
"polymer is not supported by multisegment well yet");
55 const CompletionSet& completion_set = well_ecl_->getCompletions(current_step_);
56 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
57 const Completion& completion = completion_set.get(perf);
58 const int segment_number = completion.getSegmentNumber();
59 const int segment_index = segmentNumberToIndex(segment_number);
60 segment_perforations_[segment_index].push_back(perf);
65 const Segment& segment = segmentSet()[seg];
66 const int segment_number = segment.segmentNumber();
67 const int outlet_segment_number = segment.outletSegment();
68 if (outlet_segment_number > 0) {
69 const int segment_index = segmentNumberToIndex(segment_number);
70 const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
71 segment_inlets_[outlet_segment_index].push_back(segment_index);
76 perf_depth_.resize(number_of_perforations_, 0.);
78 const double segment_depth = segmentSet()[seg].depth();
79 for (
const int perf : segment_perforations_[seg]) {
80 perf_depth_[perf] = completion_set.get(perf).getCenterDepth();
81 perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth;
88 const double segment_depth = segmentSet()[seg].depth();
89 const int outlet_segment_number = segmentSet()[seg].outletSegment();
90 const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
91 const double outlet_depth = outlet_segment.depth();
92 segment_depth_diffs_[seg] = segment_depth - outlet_depth;
100 template <
typename TypeTag>
102 MultisegmentWell<TypeTag>::
103 init(
const PhaseUsage* phase_usage_arg,
104 const std::vector<bool>* active_arg,
105 const std::vector<double>& depth_arg,
106 const double gravity_arg,
109 Base::init(phase_usage_arg, active_arg,
110 depth_arg, gravity_arg, num_cells);
122 initMatrixAndVectors(num_cells);
125 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
126 const int cell_idx = well_cells_[perf];
127 cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf];
135 template <
typename TypeTag>
137 MultisegmentWell<TypeTag>::
138 initMatrixAndVectors(
const int num_cells)
const 140 duneB_.setBuildMode( OffDiagMatWell::row_wise );
141 duneC_.setBuildMode( OffDiagMatWell::row_wise );
142 duneD_.setBuildMode( DiagMatWell::row_wise );
151 int nnz_d = numberOfSegments();
152 for (
const std::vector<int>& inlets : segment_inlets_) {
153 nnz_d += 2 * inlets.size();
155 duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d);
157 duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
158 duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
161 for (
auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) {
163 const int seg = row.index();
165 const Segment& segment = segmentSet()[seg];
166 const int outlet_segment_number = segment.outletSegment();
167 if (outlet_segment_number > 0) {
168 const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
169 row.insert(outlet_segment_index);
176 for (
const int& inlet : segment_inlets_[seg]) {
182 for (
auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) {
184 for (
const int& perf : segment_perforations_[row.index()]) {
185 const int cell_idx = well_cells_[perf];
186 row.insert(cell_idx);
191 for (
auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) {
193 for (
const int& perf : segment_perforations_[row.index()]) {
194 const int cell_idx = well_cells_[perf];
195 row.insert(cell_idx);
199 resWell_.resize( numberOfSegments() );
201 primary_variables_.resize(numberOfSegments());
202 primary_variables_evaluation_.resize(numberOfSegments());
209 template <
typename TypeTag>
211 MultisegmentWell<TypeTag>::
212 initPrimaryVariablesEvaluation()
const 214 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
215 for (
int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
216 primary_variables_evaluation_[seg][eq_idx] = 0.0;
217 primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
218 primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + numEq, 1.0);
227 template <
typename TypeTag>
229 MultisegmentWell<TypeTag>::
230 assembleWellEq(Simulator& ebosSimulator,
232 WellState& well_state,
236 const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
237 if (use_inner_iterations) {
238 iterateWellEquations(ebosSimulator, dt, well_state);
241 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
248 template <
typename TypeTag>
256 const double target = well_controls_iget_target(well_controls_, current);
257 const double* distr = well_controls_iget_distr(well_controls_, current);
258 switch (well_controls_iget_type(well_controls_, current)) {
260 well_state.bhp()[index_of_well_] = target;
261 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
262 well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_];
269 well_state.thp()[index_of_well_] = target;
306 int numPhasesWithTargetsUnderThisControl = 0;
307 for (
int phase = 0; phase < number_of_phases_; ++phase) {
308 if (distr[phase] > 0.0) {
309 numPhasesWithTargetsUnderThisControl += 1;
313 assert(numPhasesWithTargetsUnderThisControl > 0);
315 if (well_type_ == INJECTOR) {
318 assert(numPhasesWithTargetsUnderThisControl == 1);
320 for (
int phase = 0; phase < number_of_phases_; ++phase) {
321 if (distr[phase] > 0.) {
322 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target / distr[phase];
324 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = 0.;
328 initSegmentRatesWithWellRates(well_state);
330 }
else if (well_type_ == PRODUCER) {
334 double original_rates_under_phase_control = 0.0;
335 for (
int phase = 0; phase < number_of_phases_; ++phase) {
336 if (distr[phase] > 0.0) {
337 original_rates_under_phase_control += well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] * distr[phase];
341 if (original_rates_under_phase_control != 0.0 ) {
342 double scaling_factor = target / original_rates_under_phase_control;
344 for (
int phase = 0; phase < number_of_phases_; ++phase) {
345 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] *= scaling_factor;
348 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
349 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
350 well_state.segRates()[number_of_phases_ * (seg + top_segment_index) + phase] *= scaling_factor;
355 const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
356 for (
int phase = 0; phase < number_of_phases_; ++phase) {
357 if (distr[phase] > 0.0) {
358 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided / distr[phase];
361 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided;
364 initSegmentRatesWithWellRates(well_state);
371 updatePrimaryVariables(well_state);
378 template <
typename TypeTag>
383 for (
int phase = 0; phase < number_of_phases_; ++phase) {
384 const double perf_phaserate = well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] / number_of_perforations_;
385 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
386 well_state.perfPhaseRates()[number_of_phases_ * (first_perf_ + perf) + phase] = perf_phaserate;
390 const std::vector<double> perforation_rates(well_state.perfPhaseRates().begin() + number_of_phases_ * first_perf_,
391 well_state.perfPhaseRates().begin() +
392 number_of_phases_ * (first_perf_ + number_of_perforations_) );
393 std::vector<double> segment_rates;
394 WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_,
396 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
397 std::copy(segment_rates.begin(), segment_rates.end(),
398 well_state.segRates().begin() + number_of_phases_ * top_segment_index );
406 template <
typename TypeTag>
407 typename MultisegmentWell<TypeTag>::ConvergenceReport
412 assert( (
int(B_avg.size()) == numComponents()) );
415 std::vector<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
416 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
417 for (
int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
418 abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
422 std::vector<double> maximum_residual(numWellEq, 0.0);
426 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
427 for (
int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
428 if (eq_idx < numComponents()) {
429 const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
431 if (std::isnan(flux_residual)) {
432 report.nan_residual_found =
true;
433 const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
434 const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
435 report.nan_residual_wells.push_back(problem_well);
436 }
else if (flux_residual > param_.max_residual_allowed_) {
437 report.too_large_residual_found =
true;
438 const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
439 const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
440 report.nan_residual_wells.push_back(problem_well);
442 if (flux_residual > maximum_residual[eq_idx]) {
443 maximum_residual[eq_idx] = flux_residual;
449 const double pressure_residal = abs_residual[seg][eq_idx];
450 const std::string eq_name(
"Pressure");
451 if (std::isnan(pressure_residal)) {
452 report.nan_residual_found =
true;
453 const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
454 report.nan_residual_wells.push_back(problem_well);
455 }
else if (std::isinf(pressure_residal)) {
456 report.too_large_residual_found =
true;
457 const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
458 report.nan_residual_wells.push_back(problem_well);
460 if (pressure_residal > maximum_residual[eq_idx]) {
461 maximum_residual[eq_idx] = pressure_residal;
468 if ( !(report.nan_residual_found || report.too_large_residual_found) ) {
470 for (
int comp_idx = 0; comp_idx < numComponents(); ++comp_idx)
472 report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_);
475 report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_);
477 report.converged =
false;
487 template <
typename TypeTag>
490 apply(
const BVector& x, BVector& Ax)
const 492 BVectorWell Bx(duneB_.N());
498 const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx);
500 const BVectorWell invDBx = mswellhelpers::invDX(duneD_, Bx);
504 duneC_.mmtv(invDBx,Ax);
511 template <
typename TypeTag>
518 const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_);
520 const BVectorWell invDrw = mswellhelpers::invDX(duneD_, resWell_);
521 #endif // HAVE_UMFPACK 523 duneC_.mmtv(invDrw, r);
530 template <
typename TypeTag>
537 recoverSolutionWell(x, xw);
538 updateWellState(xw,
false, well_state);
545 template <
typename TypeTag>
550 std::vector<double>& )
552 OPM_THROW(std::runtime_error,
"well potential calculation for multisegment wells is not supported yet");
559 template <
typename TypeTag>
568 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
569 const std::vector<double>& segment_rates = well_state.segRates();
570 const PhaseUsage& pu = phaseUsage();
572 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
574 double total_seg_rate = 0.0;
575 const int seg_index = top_segment_index + seg;
577 primary_variables_[seg][SPres] = well_state.segPress()[seg_index];
580 for (
int p = 0; p < number_of_phases_; p++) {
581 total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p];
584 primary_variables_[seg][GTotal] = total_seg_rate;
585 if (std::abs(total_seg_rate) > 0.) {
586 if (active()[Water]) {
587 const int water_pos = pu.phase_pos[Water];
588 primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
591 const int gas_pos = pu.phase_pos[Gas];
592 primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
595 if (well_type_ == INJECTOR) {
597 const double* distr = well_controls_get_current_distr(well_controls_);
598 if (active()[Water]) {
599 if (distr[pu.phase_pos[Water]] > 0.0) {
600 primary_variables_[seg][WFrac] = 1.0;
602 primary_variables_[seg][WFrac] = 0.0;
607 if (distr[pu.phase_pos[Gas]] > 0.0) {
608 primary_variables_[seg][GFrac] = 1.0;
610 primary_variables_[seg][GFrac] = 0.0;
613 }
else if (well_type_ == PRODUCER) {
614 if (active()[Water]) {
615 primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
619 primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
630 template <
typename TypeTag>
632 MultisegmentWell<TypeTag>::
633 recoverSolutionWell(
const BVector& x, BVectorWell& xw)
const 635 BVectorWell resWell = resWell_;
637 duneB_.mmv(x, resWell);
640 xw = mswellhelpers::invDXDirect(duneD_, resWell);
642 xw = mswellhelpers::invDX(duneD_, resWell);
643 #endif // HAVE_UMFPACK 650 template <
typename TypeTag>
652 MultisegmentWell<TypeTag>::
653 solveEqAndUpdateWellState(WellState& well_state)
658 const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
660 const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
663 updateWellState(dx_well,
false, well_state);
670 template <
typename TypeTag>
672 MultisegmentWell<TypeTag>::
673 computePerfCellPressDiffs(
const Simulator& ebosSimulator)
675 for (
int perf = 0; perf < number_of_perforations_; ++perf) {
677 std::vector<double> kr(number_of_phases_, 0.0);
678 std::vector<double> density(number_of_phases_, 0.0);
680 const int cell_idx = well_cells_[perf];
681 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
682 const auto& fs = intQuants.fluidState();
686 const PhaseUsage& pu = phaseUsage();
687 if (pu.phase_used[BlackoilPhases::Aqua]) {
688 const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
689 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
690 sum_kr += kr[water_pos];
691 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
694 if (pu.phase_used[BlackoilPhases::Liquid]) {
695 const int oil_pos = pu.phase_pos[BlackoilPhases::Liquid];
696 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
697 sum_kr += kr[oil_pos];
698 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
701 if (pu.phase_used[BlackoilPhases::Vapour]) {
702 const int gas_pos = pu.phase_pos[BlackoilPhases::Vapour];
703 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
704 sum_kr += kr[gas_pos];
705 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
708 assert(sum_kr != 0.);
711 double average_density = 0.;
712 for (
int p = 0; p < number_of_phases_; ++p) {
713 average_density += kr[p] * density[p];
715 average_density /= sum_kr;
717 cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf];
725 template <
typename TypeTag>
727 MultisegmentWell<TypeTag>::
728 computeInitialComposition()
730 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
733 for (
int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) {
734 segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value();
743 template <
typename TypeTag>
745 MultisegmentWell<TypeTag>::
746 updateWellState(
const BVectorWell& dwells,
747 const bool inner_iteration,
748 WellState& well_state)
const 750 const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
752 const double relaxation_factor = (use_inner_iterations && inner_iteration) ? 0.2 : 1.0;
754 const double dFLimit = param_.dwell_fraction_max_;
755 const double max_pressure_change = param_.max_pressure_change_ms_wells_;
756 const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
758 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
759 if (active()[ Water ]) {
760 const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
761 const double dx_limited =
sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit);
762 primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
765 if (active()[ Gas ]) {
766 const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
767 const double dx_limited =
sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit);
768 primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
772 processFractions(seg);
776 const int sign = dwells[seg][SPres] > 0.? 1 : -1;
777 const double dx_limited =
sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
778 primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited;
783 primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];
788 updateWellStateFromPrimaryVariables(well_state);
795 template <
typename TypeTag>
797 MultisegmentWell<TypeTag>::
798 calculateExplicitQuantities(
const Simulator& ebosSimulator,
801 computePerfCellPressDiffs(ebosSimulator);
802 computeInitialComposition();
809 template <
typename TypeTag>
811 MultisegmentWell<TypeTag>::
814 return well_ecl_->getSegmentSet(current_step_);
821 template <
typename TypeTag>
826 return segmentSet().numberSegment();
833 template <
typename TypeTag>
838 return segmentSet().number_of_perforations_;
845 template <
typename TypeTag>
846 WellSegment::CompPressureDropEnum
847 MultisegmentWell<TypeTag>::
848 compPressureDrop()
const 850 return segmentSet().compPressureDrop();
857 template <
typename TypeTag>
858 WellSegment::MultiPhaseModelEnum
859 MultisegmentWell<TypeTag>::
860 multiphaseModel()
const 862 return segmentSet().multiPhaseModel();
869 template <
typename TypeTag>
871 MultisegmentWell<TypeTag>::
872 segmentNumberToIndex(
const int segment_number)
const 874 return segmentSet().segmentNumberToIndex(segment_number);
881 template <
typename TypeTag>
882 typename MultisegmentWell<TypeTag>::EvalWell
883 MultisegmentWell<TypeTag>::
884 volumeFraction(
const int seg,
const int comp_idx)
const 886 const PhaseUsage& pu = phaseUsage();
888 if (active()[Water] && comp_idx == pu.phase_pos[Water]) {
889 return primary_variables_evaluation_[seg][WFrac];
892 if (active()[Gas] && comp_idx == pu.phase_pos[Gas]) {
893 return primary_variables_evaluation_[seg][GFrac];
897 EvalWell oil_fraction = 1.0;
898 if (active()[Water]) {
899 oil_fraction -= primary_variables_evaluation_[seg][WFrac];
903 oil_fraction -= primary_variables_evaluation_[seg][GFrac];
915 template <
typename TypeTag>
916 typename MultisegmentWell<TypeTag>::EvalWell
917 MultisegmentWell<TypeTag>::
918 volumeFractionScaled(
const int seg,
const int comp_idx)
const 923 const double scale = scalingFactor(comp_idx);
925 return volumeFraction(seg, comp_idx) / scale;
928 return volumeFraction(seg, comp_idx);
935 template <
typename TypeTag>
936 typename MultisegmentWell<TypeTag>::EvalWell
937 MultisegmentWell<TypeTag>::
938 surfaceVolumeFraction(
const int seg,
const int comp_idx)
const 940 EvalWell sum_volume_fraction_scaled = 0.;
941 const int num_comp = numComponents();
942 for (
int idx = 0; idx < num_comp; ++idx) {
943 sum_volume_fraction_scaled += volumeFractionScaled(seg, idx);
946 assert(sum_volume_fraction_scaled.value() != 0.);
948 return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
955 template <
typename TypeTag>
957 MultisegmentWell<TypeTag>::
958 computePerfRate(
const IntensiveQuantities& int_quants,
959 const std::vector<EvalWell>& mob_perfcells,
962 const EvalWell& segment_pressure,
963 const bool& allow_cf,
964 std::vector<EvalWell>& cq_s)
const 966 const int num_comp = numComponents();
967 std::vector<EvalWell> cmix_s(num_comp, 0.0);
970 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
971 cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
974 const auto& fs = int_quants.fluidState();
976 const EvalWell pressure_cell = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
977 const EvalWell rs = extendEval(fs.Rs());
978 const EvalWell rv = extendEval(fs.Rv());
981 std::vector<EvalWell> b_perfcells(num_comp, 0.0);
983 for (
int phase = 0; phase < number_of_phases_; ++phase) {
984 const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase);
985 b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos));
989 const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
991 const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];
995 const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff);
997 const Opm::PhaseUsage& pu = phaseUsage();
1000 if ( drawdown > 0.0) {
1002 if (!allow_cf && well_type_ == INJECTOR) {
1007 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1008 const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown);
1009 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
1012 if (active()[Oil] && active()[Gas]) {
1013 const int oilpos = pu.phase_pos[Oil];
1014 const int gaspos = pu.phase_pos[Gas];
1015 const EvalWell cq_s_oil = cq_s[oilpos];
1016 const EvalWell cq_s_gas = cq_s[gaspos];
1017 cq_s[gaspos] += rs * cq_s_oil;
1018 cq_s[oilpos] += rv * cq_s_gas;
1022 if (!allow_cf && well_type_ == PRODUCER) {
1027 EvalWell total_mob = mob_perfcells[0];
1028 for (
int comp_idx = 1; comp_idx < num_comp; ++comp_idx) {
1029 total_mob += mob_perfcells[comp_idx];
1033 const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown);
1036 EvalWell volume_ratio = 0.0;
1037 if (active()[Water]) {
1038 const int watpos = pu.phase_pos[Water];
1039 volume_ratio += cmix_s[watpos] / b_perfcells[watpos];
1042 if (active()[Oil] && active()[Gas]) {
1043 const int oilpos = pu.phase_pos[Oil];
1044 const int gaspos = pu.phase_pos[Gas];
1049 const EvalWell d = 1.0 - rv * rs;
1051 if (d.value() == 0.0) {
1052 OPM_THROW(Opm::NumericalProblem,
"Zero d value obtained for well " << name() <<
" during flux calcuation" 1053 <<
" with rs " << rs <<
" and rv " << rv);
1056 const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
1057 volume_ratio += tmp_oil / b_perfcells[oilpos];
1059 const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
1060 volume_ratio += tmp_gas / b_perfcells[gaspos];
1062 if (active()[Oil]) {
1063 const int oilpos = pu.phase_pos[Oil];
1064 volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos];
1066 if (active()[Gas]) {
1067 const int gaspos = pu.phase_pos[Gas];
1068 volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos];
1072 EvalWell cqt_is = cqt_i / volume_ratio;
1073 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1074 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
1083 template <
typename TypeTag>
1084 typename MultisegmentWell<TypeTag>::EvalWell
1085 MultisegmentWell<TypeTag>::
1086 extendEval(
const Eval& in)
const 1089 out.setValue(in.value());
1090 for(
int eq_idx = 0; eq_idx < numEq;++eq_idx) {
1091 out.setDerivative(eq_idx, in.derivative(eq_idx));
1100 template <
typename TypeTag>
1102 MultisegmentWell<TypeTag>::
1103 computeSegmentFluidProperties(
const Simulator& ebosSimulator)
1112 EvalWell temperature;
1118 int pvt_region_index;
1121 const int cell_idx = well_cells_[0];
1122 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1123 const auto& fs = intQuants.fluidState();
1124 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1125 pvt_region_index = fs.pvtRegionIndex();
1128 std::vector<double> surf_dens(number_of_phases_);
1131 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1132 surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index );
1135 const int num_comp = numComponents();
1136 const Opm::PhaseUsage& pu = phaseUsage();
1137 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
1139 std::vector<EvalWell> mix_s(num_comp, 0.0);
1140 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1141 mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
1144 std::vector<EvalWell> b(num_comp, 0.0);
1146 std::vector<EvalWell> visc(number_of_phases_, 0.0);
1147 const EvalWell seg_pressure = getSegmentPressure(seg);
1148 if (pu.phase_used[BlackoilPhases::Aqua]) {
1150 const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
1152 FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1154 FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure);
1159 if (pu.phase_used[BlackoilPhases::Vapour]) {
1160 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1161 if (pu.phase_used[BlackoilPhases::Liquid]) {
1162 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1163 const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
1164 if (mix_s[oilpos] > 0.0) {
1165 if (mix_s[gaspos] > 0.0) {
1166 rv = mix_s[oilpos] / mix_s[gaspos];
1173 FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
1175 FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
1178 FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1180 FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1185 FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1187 FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1193 if (pu.phase_used[BlackoilPhases::Liquid]) {
1194 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1195 if (pu.phase_used[BlackoilPhases::Liquid]) {
1196 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1197 const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
1198 if (mix_s[gaspos] > 0.0) {
1199 if (mix_s[oilpos] > 0.0) {
1200 rs = mix_s[gaspos] / mix_s[oilpos];
1207 FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
1209 FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
1212 FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1214 FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1219 FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
1221 FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
1225 std::vector<EvalWell> mix(mix_s);
1226 if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
1227 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
1228 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
1230 mix[gaspos] = (mix_s[gaspos] - mix_s[oilpos] * rs) / (1. - rs * rv);
1233 mix[oilpos] = (mix_s[oilpos] - mix_s[gaspos] * rv) / (1. - rs * rv);
1237 EvalWell volrat(0.0);
1238 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1239 volrat += mix[comp_idx] / b[comp_idx];
1242 segment_viscosities_[seg] = 0.;
1244 for (
int p = 0; p < number_of_phases_; ++p) {
1245 const EvalWell phase_fraction = mix[p] / b[p] / volrat;
1246 segment_viscosities_[seg] += visc[p] * phase_fraction;
1249 EvalWell density(0.0);
1250 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1251 density += surf_dens[comp_idx] * mix_s[comp_idx];
1253 segment_densities_[seg] = density / volrat;
1256 segment_mass_rates_[seg] = 0.;
1257 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1258 const EvalWell rate = getSegmentRate(seg, phase);
1259 segment_mass_rates_[seg] += rate * surf_dens[phase];
1268 template <
typename TypeTag>
1269 typename MultisegmentWell<TypeTag>::EvalWell
1270 MultisegmentWell<TypeTag>::
1271 getSegmentPressure(
const int seg)
const 1273 return primary_variables_evaluation_[seg][SPres];
1280 template <
typename TypeTag>
1281 typename MultisegmentWell<TypeTag>::EvalWell
1282 MultisegmentWell<TypeTag>::
1283 getSegmentRate(
const int seg,
1284 const int comp_idx)
const 1286 return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx);
1293 template <
typename TypeTag>
1294 typename MultisegmentWell<TypeTag>::EvalWell
1295 MultisegmentWell<TypeTag>::
1296 getSegmentGTotal(
const int seg)
const 1298 return primary_variables_evaluation_[seg][GTotal];
1305 template <
typename TypeTag>
1307 MultisegmentWell<TypeTag>::
1308 getMobility(
const Simulator& ebosSimulator,
1310 std::vector<EvalWell>& mob)
const 1313 const int np = number_of_phases_;
1314 const int cell_idx = well_cells_[perf];
1315 assert (
int(mob.size()) == numComponents());
1316 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1317 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1321 const int satid = saturation_table_number_[perf] - 1;
1322 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1323 if( satid == satid_elem ) {
1325 for (
int phase = 0; phase < np; ++phase) {
1326 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
1327 mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
1334 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1335 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
1336 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1339 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1342 for (
int phase = 0; phase < np; ++phase) {
1343 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
1344 mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
1354 template <
typename TypeTag>
1356 MultisegmentWell<TypeTag>::
1357 assembleControlEq()
const 1359 EvalWell control_eq(0.0);
1361 switch (well_controls_get_current_type(well_controls_)) {
1364 OPM_THROW(std::runtime_error,
"Not handling THP control for Multisegment wells for now");
1368 const double target_bhp = well_controls_get_current_target(well_controls_);
1369 control_eq = getSegmentPressure(0) - target_bhp;
1375 int number_phases_under_control = 0;
1376 const double* distr = well_controls_get_current_distr(well_controls_);
1377 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1378 if (distr[phase] > 0.0) {
1379 ++number_phases_under_control;
1382 assert(number_phases_under_control > 0);
1384 const std::vector<double> g = {1.0, 1.0, 0.01};
1385 const double target_rate = well_controls_get_current_target(well_controls_);
1388 if (number_phases_under_control == 1) {
1389 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1390 if (distr[phase] > 0.) {
1391 control_eq = getSegmentGTotal(0) * volumeFraction(0, phase) - g[phase] * target_rate;
1396 EvalWell rate_for_control(0.0);
1397 const EvalWell G_total = getSegmentGTotal(0);
1398 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1399 if (distr[phase] > 0.) {
1400 rate_for_control += G_total * volumeFractionScaled(0, phase);
1404 control_eq = rate_for_control - target_rate;
1408 case RESERVOIR_RATE:
1410 EvalWell rate_for_control(0.0);
1411 const double* distr = well_controls_get_current_distr(well_controls_);
1412 for (
int phase = 0; phase < number_of_phases_; ++phase) {
1413 if (distr[phase] > 0.) {
1414 rate_for_control += getSegmentGTotal(0) * volumeFraction(0, phase);
1417 const double target_rate = well_controls_get_current_target(well_controls_);
1418 control_eq = rate_for_control - target_rate;
1422 OPM_THROW(std::runtime_error,
"Unknown well control control types for well " << name());
1428 resWell_[0][SPres] = control_eq.value();
1429 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1430 duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq);
1438 template <
typename TypeTag>
1440 MultisegmentWell<TypeTag>::
1441 assemblePressureEq(
const int seg)
const 1446 EvalWell pressure_equation = getSegmentPressure(seg);
1450 pressure_equation -= getHydroPressureLoss(seg);
1452 if (frictionalPressureLossConsidered()) {
1453 pressure_equation -= getFrictionPressureLoss(seg);
1456 resWell_[seg][SPres] = pressure_equation.value();
1457 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1458 duneD_[seg][seg][SPres][pv_idx] = pressure_equation.derivative(pv_idx + numEq);
1462 const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
1463 const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
1465 resWell_[seg][SPres] -= outlet_pressure.value();
1466 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1467 duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
1470 if (accelerationalPressureLossConsidered()) {
1471 handleAccelerationPressureLoss(seg);
1479 template <
typename TypeTag>
1480 typename MultisegmentWell<TypeTag>::EvalWell
1481 MultisegmentWell<TypeTag>::
1482 getHydroPressureLoss(
const int seg)
const 1484 return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg];
1491 template <
typename TypeTag>
1492 typename MultisegmentWell<TypeTag>::EvalWell
1493 MultisegmentWell<TypeTag>::
1494 getFrictionPressureLoss(
const int seg)
const 1496 const EvalWell mass_rate = segment_mass_rates_[seg];
1497 const EvalWell density = segment_densities_[seg];
1498 const EvalWell visc = segment_viscosities_[seg];
1499 const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
1500 const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength();
1501 assert(length > 0.);
1502 const double roughness = segmentSet()[seg].roughness();
1503 const double area = segmentSet()[seg].crossArea();
1504 const double diameter = segmentSet()[seg].internalDiameter();
1506 const double sign = mass_rate < 0. ? 1.0 : - 1.0;
1508 return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
1515 template <
typename TypeTag>
1517 MultisegmentWell<TypeTag>::
1518 handleAccelerationPressureLoss(
const int seg)
const 1522 const double area = segmentSet()[seg].crossArea();
1523 const EvalWell mass_rate = segment_mass_rates_[seg];
1524 const EvalWell density = segment_densities_[seg];
1525 const EvalWell out_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
1527 resWell_[seg][SPres] -= out_velocity_head.value();
1528 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1529 duneD_[seg][seg][SPres][pv_idx] -= out_velocity_head.derivative(pv_idx + numEq);
1533 double max_area = area;
1534 for (
const int inlet : segment_inlets_[seg]) {
1535 const double inlet_area = segmentSet()[inlet].crossArea();
1536 if (inlet_area > max_area) {
1537 max_area = inlet_area;
1542 for (
const int inlet : segment_inlets_[seg]) {
1543 const EvalWell density = segment_densities_[inlet];
1544 const EvalWell mass_rate = segment_mass_rates_[inlet];
1545 const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
1546 resWell_[seg][SPres] += inlet_velocity_head.value();
1547 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1548 duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq);
1557 template <
typename TypeTag>
1559 MultisegmentWell<TypeTag>::
1560 processFractions(
const int seg)
const 1562 const PhaseUsage& pu = phaseUsage();
1564 std::vector<double> fractions(number_of_phases_, 0.0);
1566 assert( active()[Oil] );
1567 const int oil_pos = pu.phase_pos[Oil];
1568 fractions[oil_pos] = 1.0;
1570 if ( active()[Water] ) {
1571 const int water_pos = pu.phase_pos[Water];
1572 fractions[water_pos] = primary_variables_[seg][WFrac];
1573 fractions[oil_pos] -= fractions[water_pos];
1576 if ( active()[Gas] ) {
1577 const int gas_pos = pu.phase_pos[Gas];
1578 fractions[gas_pos] = primary_variables_[seg][GFrac];
1579 fractions[oil_pos] -= fractions[gas_pos];
1582 if (active()[Water]) {
1583 const int water_pos = pu.phase_pos[Water];
1584 if (fractions[water_pos] < 0.0) {
1585 if ( active()[Gas] ) {
1586 fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
1588 fractions[oil_pos] /= (1.0 - fractions[water_pos]);
1589 fractions[water_pos] = 0.0;
1593 if (active()[Gas]) {
1594 const int gas_pos = pu.phase_pos[Gas];
1595 if (fractions[gas_pos] < 0.0) {
1596 if ( active()[Water] ) {
1597 fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
1599 fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
1600 fractions[gas_pos] = 0.0;
1604 if (fractions[oil_pos] < 0.0) {
1605 if ( active()[Water] ) {
1606 fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
1608 if ( active()[Gas] ) {
1609 fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
1611 fractions[oil_pos] = 0.0;
1614 if ( active()[Water] ) {
1615 primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
1618 if ( active()[Gas] ) {
1619 primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
1627 template <
typename TypeTag>
1629 MultisegmentWell<TypeTag>::
1630 updateWellStateFromPrimaryVariables(WellState& well_state)
const 1632 const PhaseUsage& pu = phaseUsage();
1633 assert( active()[Oil] );
1634 const int oil_pos = pu.phase_pos[Oil];
1636 for (
int seg = 0; seg < numberOfSegments(); ++seg) {
1637 std::vector<double> fractions(number_of_phases_, 0.0);
1638 fractions[oil_pos] = 1.0;
1640 if ( active()[Water] ) {
1641 const int water_pos = pu.phase_pos[Water];
1642 fractions[water_pos] = primary_variables_[seg][WFrac];
1643 fractions[oil_pos] -= fractions[water_pos];
1646 if ( active()[Gas] ) {
1647 const int gas_pos = pu.phase_pos[Gas];
1648 fractions[gas_pos] = primary_variables_[seg][GFrac];
1649 fractions[oil_pos] -= fractions[gas_pos];
1653 for (
int p = 0; p < number_of_phases_; ++p) {
1654 const double scale = scalingFactor(p);
1657 fractions[p] /= scale;
1665 const double g_total = primary_variables_[seg][GTotal];
1666 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
1667 for (
int p = 0; p < number_of_phases_; ++p) {
1668 const double phase_rate = g_total * fractions[p];
1669 well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate;
1671 well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate;
1676 well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres];
1678 well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index];
1687 template <
typename TypeTag>
1689 MultisegmentWell<TypeTag>::
1690 scalingFactor(
const int comp_idx)
const 1692 const double* distr = well_controls_get_current_distr(well_controls_);
1694 if (well_controls_get_current_type(well_controls_) == RESERVOIR_RATE) {
1697 return distr[comp_idx];
1700 const PhaseUsage& pu = phaseUsage();
1702 if (active()[Water] && pu.phase_pos[Water] == comp_idx)
1704 if (active()[Oil] && pu.phase_pos[Oil] == comp_idx)
1706 if (active()[Gas] && pu.phase_pos[Gas] == comp_idx)
1720 template <
typename TypeTag>
1722 MultisegmentWell<TypeTag>::
1723 frictionalPressureLossConsidered()
const 1726 return (segmentSet().compPressureDrop() != WellSegment::H__);
1733 template <
typename TypeTag>
1735 MultisegmentWell<TypeTag>::
1736 accelerationalPressureLossConsidered()
const 1738 return (segmentSet().compPressureDrop() == WellSegment::HFA);
1745 template<
typename TypeTag>
1747 MultisegmentWell<TypeTag>::
1748 iterateWellEquations(Simulator& ebosSimulator,
1750 WellState& well_state)
1757 const int max_iter_number = param_.max_inner_iter_ms_wells_;
1759 for (; it < max_iter_number; ++it) {
1761 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state,
true);
1764 const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
1766 const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
1774 const std::vector<double> B {0.5, 0.5, 0.005};
1776 const ConvergenceReport report = getWellConvergence(B);
1778 if (report.converged) {
1782 updateWellState(dx_well,
true, well_state);
1784 initPrimaryVariablesEvaluation();
1793 template<
typename TypeTag>
1795 MultisegmentWell<TypeTag>::
1796 assembleWellEqWithoutIteration(Simulator& ebosSimulator,
1798 WellState& well_state,
1802 computeSegmentFluidProperties(ebosSimulator);
1818 auto& ebosJac = ebosSimulator.model().linearizer().matrix();
1819 auto& ebosResid = ebosSimulator.model().linearizer().residual();
1821 const bool allow_cf = getAllowCrossFlow();
1823 const int nseg = numberOfSegments();
1824 const int num_comp = numComponents();
1826 for (
int seg = 0; seg < nseg; ++seg) {
1830 const double volume = segmentSet()[seg].volume();
1832 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1833 EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx])
1834 + getSegmentRate(seg, comp_idx);
1836 resWell_[seg][comp_idx] += accumulation_term.value();
1837 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1838 duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
1845 for (
const int inlet : segment_inlets_[seg]) {
1846 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1847 const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx);
1848 resWell_[seg][comp_idx] -= inlet_rate.value();
1849 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1850 duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq);
1857 const EvalWell seg_pressure = getSegmentPressure(seg);
1858 for (
const int perf : segment_perforations_[seg]) {
1859 const int cell_idx = well_cells_[perf];
1860 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
1861 std::vector<EvalWell> mob(num_comp, 0.0);
1862 getMobility(ebosSimulator, perf, mob);
1863 std::vector<EvalWell> cq_s(num_comp, 0.0);
1864 computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s);
1866 for (
int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
1868 const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
1875 ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
1879 resWell_[seg][comp_idx] -= cq_s_effective.value();
1882 for (
int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1885 duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq);
1888 duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
1891 for (
int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
1894 ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx);
1895 duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
1906 assembleControlEq();
1908 assemblePressureEq(seg);
virtual void recoverWellSolutionAndUpdateWellState(const BVector &x, WellState &well_state) const
using the solution x to recover the solution xw for wells and applying xw to update Well State ...
Definition: MultisegmentWell_impl.hpp:533
virtual ConvergenceReport getWellConvergence(const std::vector< double > &B_avg) const
check whether the well equations get converged for this well
Definition: MultisegmentWell_impl.hpp:409
Definition: AutoDiff.hpp:297
a struct to collect information about the convergence checking
Definition: WellInterface.hpp:117
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
int numberOfSegments() const
number of segments for this well int number_of_segments_;
Definition: MultisegmentWell_impl.hpp:824
virtual void apply(const BVector &x, BVector &Ax) const
Ax = Ax - C D^-1 B x.
Definition: MultisegmentWell_impl.hpp:490
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
virtual void computeWellPotentials(const Simulator &ebosSimulator, const WellState &well_state, std::vector< double > &well_potentials)
computing the well potentials for group control
Definition: MultisegmentWell_impl.hpp:548
virtual void updateWellStateWithTarget(const int current, WellState &well_state) const
updating the well state based the control mode specified with current
Definition: MultisegmentWell_impl.hpp:251
Definition: MultisegmentWell.hpp:32
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719