00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <opm/autodiff/MSWellHelpers.hpp>
00023
00024 namespace Opm
00025 {
00026
00027
00028 template <typename TypeTag>
00029 MultisegmentWell<TypeTag>::
00030 MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
00031 : Base(well, time_step, wells, param)
00032 , segment_perforations_(numberOfSegments())
00033 , segment_inlets_(numberOfSegments())
00034 , cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
00035 , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0)
00036 , perforation_segment_depth_diffs_(number_of_perforations_, 0.0)
00037 , segment_comp_initial_(numberOfSegments(), std::vector<double>(numComponents(), 0.0))
00038 , segment_densities_(numberOfSegments(), 0.0)
00039 , segment_viscosities_(numberOfSegments(), 0.0)
00040 , segment_mass_rates_(numberOfSegments(), 0.0)
00041 , segment_depth_diffs_(numberOfSegments(), 0.0)
00042 {
00043
00044 if (has_solvent) {
00045 OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
00046 }
00047
00048 if (has_polymer) {
00049 OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
00050 }
00051
00052
00053
00054
00055 const CompletionSet& completion_set = well_ecl_->getCompletions(current_step_);
00056 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00057 const Completion& completion = completion_set.get(perf);
00058 const int segment_number = completion.getSegmentNumber();
00059 const int segment_index = segmentNumberToIndex(segment_number);
00060 segment_perforations_[segment_index].push_back(perf);
00061 }
00062
00063
00064 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00065 const Segment& segment = segmentSet()[seg];
00066 const int segment_number = segment.segmentNumber();
00067 const int outlet_segment_number = segment.outletSegment();
00068 if (outlet_segment_number > 0) {
00069 const int segment_index = segmentNumberToIndex(segment_number);
00070 const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
00071 segment_inlets_[outlet_segment_index].push_back(segment_index);
00072 }
00073 }
00074
00075
00076 perf_depth_.resize(number_of_perforations_, 0.);
00077 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00078 const double segment_depth = segmentSet()[seg].depth();
00079 for (const int perf : segment_perforations_[seg]) {
00080 perf_depth_[perf] = completion_set.get(perf).getCenterDepth();
00081 perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth;
00082 }
00083 }
00084
00085
00086
00087 for (int seg = 1; seg < numberOfSegments(); ++seg) {
00088 const double segment_depth = segmentSet()[seg].depth();
00089 const int outlet_segment_number = segmentSet()[seg].outletSegment();
00090 const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
00091 const double outlet_depth = outlet_segment.depth();
00092 segment_depth_diffs_[seg] = segment_depth - outlet_depth;
00093 }
00094 }
00095
00096
00097
00098
00099
00100 template <typename TypeTag>
00101 void
00102 MultisegmentWell<TypeTag>::
00103 init(const PhaseUsage* phase_usage_arg,
00104 const std::vector<bool>* active_arg,
00105 const std::vector<double>& depth_arg,
00106 const double gravity_arg,
00107 const int num_cells)
00108 {
00109 Base::init(phase_usage_arg, active_arg,
00110 depth_arg, gravity_arg, num_cells);
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 initMatrixAndVectors(num_cells);
00123
00124
00125 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00126 const int cell_idx = well_cells_[perf];
00127 cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf];
00128 }
00129 }
00130
00131
00132
00133
00134
00135 template <typename TypeTag>
00136 void
00137 MultisegmentWell<TypeTag>::
00138 initMatrixAndVectors(const int num_cells) const
00139 {
00140 duneB_.setBuildMode( OffDiagMatWell::row_wise );
00141 duneC_.setBuildMode( OffDiagMatWell::row_wise );
00142 duneD_.setBuildMode( DiagMatWell::row_wise );
00143
00144
00145
00146
00147
00148
00149
00150 {
00151 int nnz_d = numberOfSegments();
00152 for (const std::vector<int>& inlets : segment_inlets_) {
00153 nnz_d += 2 * inlets.size();
00154 }
00155 duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d);
00156 }
00157 duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
00158 duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
00159
00160
00161 for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) {
00162
00163 const int seg = row.index();
00164
00165 const Segment& segment = segmentSet()[seg];
00166 const int outlet_segment_number = segment.outletSegment();
00167 if (outlet_segment_number > 0) {
00168 const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
00169 row.insert(outlet_segment_index);
00170 }
00171
00172
00173 row.insert(seg);
00174
00175
00176 for (const int& inlet : segment_inlets_[seg]) {
00177 row.insert(inlet);
00178 }
00179 }
00180
00181
00182 for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) {
00183
00184 for (const int& perf : segment_perforations_[row.index()]) {
00185 const int cell_idx = well_cells_[perf];
00186 row.insert(cell_idx);
00187 }
00188 }
00189
00190
00191 for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) {
00192
00193 for (const int& perf : segment_perforations_[row.index()]) {
00194 const int cell_idx = well_cells_[perf];
00195 row.insert(cell_idx);
00196 }
00197 }
00198
00199 resWell_.resize( numberOfSegments() );
00200
00201 primary_variables_.resize(numberOfSegments());
00202 primary_variables_evaluation_.resize(numberOfSegments());
00203 }
00204
00205
00206
00207
00208
00209 template <typename TypeTag>
00210 void
00211 MultisegmentWell<TypeTag>::
00212 initPrimaryVariablesEvaluation() const
00213 {
00214 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00215 for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
00216 primary_variables_evaluation_[seg][eq_idx] = 0.0;
00217 primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
00218 primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + numEq, 1.0);
00219 }
00220 }
00221 }
00222
00223
00224
00225
00226
00227 template <typename TypeTag>
00228 void
00229 MultisegmentWell<TypeTag>::
00230 assembleWellEq(Simulator& ebosSimulator,
00231 const double dt,
00232 WellState& well_state,
00233 bool only_wells)
00234 {
00235
00236 const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
00237 if (use_inner_iterations) {
00238 iterateWellEquations(ebosSimulator, dt, well_state);
00239 }
00240
00241 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
00242 }
00243
00244
00245
00246
00247
00248 template <typename TypeTag>
00249 void
00250 MultisegmentWell<TypeTag>::
00251 updateWellStateWithTarget(const int current,
00252 WellState& well_state) const
00253 {
00254
00255
00256 const double target = well_controls_iget_target(well_controls_, current);
00257 const double* distr = well_controls_iget_distr(well_controls_, current);
00258 switch (well_controls_iget_type(well_controls_, current)) {
00259 case BHP: {
00260 well_state.bhp()[index_of_well_] = target;
00261 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
00262 well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_];
00263
00264
00265 break;
00266 }
00267
00268 case THP: {
00269 well_state.thp()[index_of_well_] = target;
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 break;
00293 }
00294
00295 case RESERVOIR_RATE:
00296 case SURFACE_RATE:
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 int numPhasesWithTargetsUnderThisControl = 0;
00307 for (int phase = 0; phase < number_of_phases_; ++phase) {
00308 if (distr[phase] > 0.0) {
00309 numPhasesWithTargetsUnderThisControl += 1;
00310 }
00311 }
00312
00313 assert(numPhasesWithTargetsUnderThisControl > 0);
00314
00315 if (well_type_ == INJECTOR) {
00316
00317
00318 assert(numPhasesWithTargetsUnderThisControl == 1);
00319
00320 for (int phase = 0; phase < number_of_phases_; ++phase) {
00321 if (distr[phase] > 0.) {
00322 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target / distr[phase];
00323 } else {
00324 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = 0.;
00325 }
00326 }
00327
00328 initSegmentRatesWithWellRates(well_state);
00329
00330 } else if (well_type_ == PRODUCER) {
00331
00332
00333
00334 double original_rates_under_phase_control = 0.0;
00335 for (int phase = 0; phase < number_of_phases_; ++phase) {
00336 if (distr[phase] > 0.0) {
00337 original_rates_under_phase_control += well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] * distr[phase];
00338 }
00339 }
00340
00341 if (original_rates_under_phase_control != 0.0 ) {
00342 double scaling_factor = target / original_rates_under_phase_control;
00343
00344 for (int phase = 0; phase < number_of_phases_; ++phase) {
00345 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] *= scaling_factor;
00346
00347
00348 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
00349 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00350 well_state.segRates()[number_of_phases_ * (seg + top_segment_index) + phase] *= scaling_factor;
00351 }
00352 }
00353 } else {
00354
00355 const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
00356 for (int phase = 0; phase < number_of_phases_; ++phase) {
00357 if (distr[phase] > 0.0) {
00358 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided / distr[phase];
00359 } else {
00360
00361 well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided;
00362 }
00363 }
00364 initSegmentRatesWithWellRates(well_state);
00365 }
00366 }
00367
00368 break;
00369 }
00370
00371 updatePrimaryVariables(well_state);
00372 }
00373
00374
00375
00376
00377
00378 template <typename TypeTag>
00379 void
00380 MultisegmentWell<TypeTag>::
00381 initSegmentRatesWithWellRates(WellState& well_state) const
00382 {
00383 for (int phase = 0; phase < number_of_phases_; ++phase) {
00384 const double perf_phaserate = well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] / number_of_perforations_;
00385 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00386 well_state.perfPhaseRates()[number_of_phases_ * (first_perf_ + perf) + phase] = perf_phaserate;
00387 }
00388 }
00389
00390 const std::vector<double> perforation_rates(well_state.perfPhaseRates().begin() + number_of_phases_ * first_perf_,
00391 well_state.perfPhaseRates().begin() +
00392 number_of_phases_ * (first_perf_ + number_of_perforations_) );
00393 std::vector<double> segment_rates;
00394 WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_,
00395 0, segment_rates);
00396 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
00397 std::copy(segment_rates.begin(), segment_rates.end(),
00398 well_state.segRates().begin() + number_of_phases_ * top_segment_index );
00399
00400 }
00401
00402
00403
00404
00405
00406 template <typename TypeTag>
00407 typename MultisegmentWell<TypeTag>::ConvergenceReport
00408 MultisegmentWell<TypeTag>::
00409 getWellConvergence(const std::vector<double>& B_avg) const
00410 {
00411
00412 assert( (int(B_avg.size()) == numComponents()) );
00413
00414
00415 std::vector<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
00416 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00417 for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
00418 abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
00419 }
00420 }
00421
00422 std::vector<double> maximum_residual(numWellEq, 0.0);
00423
00424 ConvergenceReport report;
00425
00426 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00427 for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
00428 if (eq_idx < numComponents()) {
00429 const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
00430
00431 if (std::isnan(flux_residual)) {
00432 report.nan_residual_found = true;
00433 const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
00434 const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
00435 report.nan_residual_wells.push_back(problem_well);
00436 } else if (flux_residual > param_.max_residual_allowed_) {
00437 report.too_large_residual_found = true;
00438 const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
00439 const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
00440 report.nan_residual_wells.push_back(problem_well);
00441 } else {
00442 if (flux_residual > maximum_residual[eq_idx]) {
00443 maximum_residual[eq_idx] = flux_residual;
00444 }
00445 }
00446 } else {
00447
00448
00449 const double pressure_residal = abs_residual[seg][eq_idx];
00450 const std::string eq_name("Pressure");
00451 if (std::isnan(pressure_residal)) {
00452 report.nan_residual_found = true;
00453 const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
00454 report.nan_residual_wells.push_back(problem_well);
00455 } else if (std::isinf(pressure_residal)) {
00456 report.too_large_residual_found = true;
00457 const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
00458 report.nan_residual_wells.push_back(problem_well);
00459 } else {
00460 if (pressure_residal > maximum_residual[eq_idx]) {
00461 maximum_residual[eq_idx] = pressure_residal;
00462 }
00463 }
00464 }
00465 }
00466 }
00467
00468 if ( !(report.nan_residual_found || report.too_large_residual_found) ) {
00469
00470 for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx)
00471 {
00472 report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_);
00473 }
00474
00475 report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_);
00476 } else {
00477 report.converged = false;
00478 }
00479
00480 return report;
00481 }
00482
00483
00484
00485
00486
00487 template <typename TypeTag>
00488 void
00489 MultisegmentWell<TypeTag>::
00490 apply(const BVector& x, BVector& Ax) const
00491 {
00492 BVectorWell Bx(duneB_.N());
00493
00494 duneB_.mv(x, Bx);
00495
00496
00497 #if HAVE_UMFPACK
00498 const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx);
00499 #else
00500 const BVectorWell invDBx = mswellhelpers::invDX(duneD_, Bx);
00501 #endif
00502
00503
00504 duneC_.mmtv(invDBx,Ax);
00505 }
00506
00507
00508
00509
00510
00511 template <typename TypeTag>
00512 void
00513 MultisegmentWell<TypeTag>::
00514 apply(BVector& r) const
00515 {
00516
00517 #if HAVE_UMFPACK
00518 const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_);
00519 #else
00520 const BVectorWell invDrw = mswellhelpers::invDX(duneD_, resWell_);
00521 #endif // HAVE_UMFPACK
00522
00523 duneC_.mmtv(invDrw, r);
00524 }
00525
00526
00527
00528
00529
00530 template <typename TypeTag>
00531 void
00532 MultisegmentWell<TypeTag>::
00533 recoverWellSolutionAndUpdateWellState(const BVector& x,
00534 WellState& well_state) const
00535 {
00536 BVectorWell xw(1);
00537 recoverSolutionWell(x, xw);
00538 updateWellState(xw, false, well_state);
00539 }
00540
00541
00542
00543
00544
00545 template <typename TypeTag>
00546 void
00547 MultisegmentWell<TypeTag>::
00548 computeWellPotentials(const Simulator& ,
00549 const WellState& ,
00550 std::vector<double>& )
00551 {
00552 OPM_THROW(std::runtime_error, "well potential calculation for multisegment wells is not supported yet");
00553 }
00554
00555
00556
00557
00558
00559 template <typename TypeTag>
00560 void
00561 MultisegmentWell<TypeTag>::
00562 updatePrimaryVariables(const WellState& well_state) const
00563 {
00564
00565
00566
00567
00568 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
00569 const std::vector<double>& segment_rates = well_state.segRates();
00570 const PhaseUsage& pu = phaseUsage();
00571
00572 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00573
00574 double total_seg_rate = 0.0;
00575 const int seg_index = top_segment_index + seg;
00576
00577 primary_variables_[seg][SPres] = well_state.segPress()[seg_index];
00578
00579
00580 for (int p = 0; p < number_of_phases_; p++) {
00581 total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p];
00582 }
00583
00584 primary_variables_[seg][GTotal] = total_seg_rate;
00585 if (std::abs(total_seg_rate) > 0.) {
00586 if (active()[Water]) {
00587 const int water_pos = pu.phase_pos[Water];
00588 primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
00589 }
00590 if (active()[Gas]) {
00591 const int gas_pos = pu.phase_pos[Gas];
00592 primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
00593 }
00594 } else {
00595 if (well_type_ == INJECTOR) {
00596
00597 const double* distr = well_controls_get_current_distr(well_controls_);
00598 if (active()[Water]) {
00599 if (distr[pu.phase_pos[Water]] > 0.0) {
00600 primary_variables_[seg][WFrac] = 1.0;
00601 } else {
00602 primary_variables_[seg][WFrac] = 0.0;
00603 }
00604 }
00605
00606 if (active()[Gas]) {
00607 if (distr[pu.phase_pos[Gas]] > 0.0) {
00608 primary_variables_[seg][GFrac] = 1.0;
00609 } else {
00610 primary_variables_[seg][GFrac] = 0.0;
00611 }
00612 }
00613 } else if (well_type_ == PRODUCER) {
00614 if (active()[Water]) {
00615 primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
00616 }
00617
00618 if (active()[Gas]) {
00619 primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
00620 }
00621 }
00622 }
00623 }
00624 }
00625
00626
00627
00628
00629
00630 template <typename TypeTag>
00631 void
00632 MultisegmentWell<TypeTag>::
00633 recoverSolutionWell(const BVector& x, BVectorWell& xw) const
00634 {
00635 BVectorWell resWell = resWell_;
00636
00637 duneB_.mmv(x, resWell);
00638
00639 #if HAVE_UMFPACK
00640 xw = mswellhelpers::invDXDirect(duneD_, resWell);
00641 #else
00642 xw = mswellhelpers::invDX(duneD_, resWell);
00643 #endif // HAVE_UMFPACK
00644 }
00645
00646
00647
00648
00649
00650 template <typename TypeTag>
00651 void
00652 MultisegmentWell<TypeTag>::
00653 solveEqAndUpdateWellState(WellState& well_state)
00654 {
00655 #if HAVE_UMFPACK
00656
00657
00658 const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
00659 #else
00660 const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
00661 #endif
00662
00663 updateWellState(dx_well, false, well_state);
00664 }
00665
00666
00667
00668
00669
00670 template <typename TypeTag>
00671 void
00672 MultisegmentWell<TypeTag>::
00673 computePerfCellPressDiffs(const Simulator& ebosSimulator)
00674 {
00675 for (int perf = 0; perf < number_of_perforations_; ++perf) {
00676
00677 std::vector<double> kr(number_of_phases_, 0.0);
00678 std::vector<double> density(number_of_phases_, 0.0);
00679
00680 const int cell_idx = well_cells_[perf];
00681 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
00682 const auto& fs = intQuants.fluidState();
00683
00684 double sum_kr = 0.;
00685
00686 const PhaseUsage& pu = phaseUsage();
00687 if (pu.phase_used[BlackoilPhases::Aqua]) {
00688 const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
00689 kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
00690 sum_kr += kr[water_pos];
00691 density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
00692 }
00693
00694 if (pu.phase_used[BlackoilPhases::Liquid]) {
00695 const int oil_pos = pu.phase_pos[BlackoilPhases::Liquid];
00696 kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
00697 sum_kr += kr[oil_pos];
00698 density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
00699 }
00700
00701 if (pu.phase_used[BlackoilPhases::Vapour]) {
00702 const int gas_pos = pu.phase_pos[BlackoilPhases::Vapour];
00703 kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
00704 sum_kr += kr[gas_pos];
00705 density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
00706 }
00707
00708 assert(sum_kr != 0.);
00709
00710
00711 double average_density = 0.;
00712 for (int p = 0; p < number_of_phases_; ++p) {
00713 average_density += kr[p] * density[p];
00714 }
00715 average_density /= sum_kr;
00716
00717 cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf];
00718 }
00719 }
00720
00721
00722
00723
00724
00725 template <typename TypeTag>
00726 void
00727 MultisegmentWell<TypeTag>::
00728 computeInitialComposition()
00729 {
00730 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00731
00732
00733 for (int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) {
00734 segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value();
00735 }
00736 }
00737 }
00738
00739
00740
00741
00742
00743 template <typename TypeTag>
00744 void
00745 MultisegmentWell<TypeTag>::
00746 updateWellState(const BVectorWell& dwells,
00747 const bool inner_iteration,
00748 WellState& well_state) const
00749 {
00750 const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
00751
00752 const double relaxation_factor = (use_inner_iterations && inner_iteration) ? 0.2 : 1.0;
00753
00754 const double dFLimit = param_.dwell_fraction_max_;
00755 const double max_pressure_change = param_.max_pressure_change_ms_wells_;
00756 const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
00757
00758 for (int seg = 0; seg < numberOfSegments(); ++seg) {
00759 if (active()[ Water ]) {
00760 const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
00761 const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit);
00762 primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
00763 }
00764
00765 if (active()[ Gas ]) {
00766 const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
00767 const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit);
00768 primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
00769 }
00770
00771
00772 processFractions(seg);
00773
00774
00775 {
00776 const int sign = dwells[seg][SPres] > 0.? 1 : -1;
00777 const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
00778 primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited;
00779 }
00780
00781
00782 {
00783 primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];
00784 }
00785
00786 }
00787
00788 updateWellStateFromPrimaryVariables(well_state);
00789 }
00790
00791
00792
00793
00794
00795 template <typename TypeTag>
00796 void
00797 MultisegmentWell<TypeTag>::
00798 calculateExplicitQuantities(const Simulator& ebosSimulator,
00799 const WellState& )
00800 {
00801 computePerfCellPressDiffs(ebosSimulator);
00802 computeInitialComposition();
00803 }
00804
00805
00806
00807
00808
00809 template <typename TypeTag>
00810 const SegmentSet&
00811 MultisegmentWell<TypeTag>::
00812 segmentSet() const
00813 {
00814 return well_ecl_->getSegmentSet(current_step_);
00815 }
00816
00817
00818
00819
00820
00821 template <typename TypeTag>
00822 int
00823 MultisegmentWell<TypeTag>::
00824 numberOfSegments() const
00825 {
00826 return segmentSet().numberSegment();
00827 }
00828
00829
00830
00831
00832
00833 template <typename TypeTag>
00834 int
00835 MultisegmentWell<TypeTag>::
00836 numberOfPerforations() const
00837 {
00838 return segmentSet().number_of_perforations_;
00839 }
00840
00841
00842
00843
00844
00845 template <typename TypeTag>
00846 WellSegment::CompPressureDropEnum
00847 MultisegmentWell<TypeTag>::
00848 compPressureDrop() const
00849 {
00850 return segmentSet().compPressureDrop();
00851 }
00852
00853
00854
00855
00856
00857 template <typename TypeTag>
00858 WellSegment::MultiPhaseModelEnum
00859 MultisegmentWell<TypeTag>::
00860 multiphaseModel() const
00861 {
00862 return segmentSet().multiPhaseModel();
00863 }
00864
00865
00866
00867
00868
00869 template <typename TypeTag>
00870 int
00871 MultisegmentWell<TypeTag>::
00872 segmentNumberToIndex(const int segment_number) const
00873 {
00874 return segmentSet().segmentNumberToIndex(segment_number);
00875 }
00876
00877
00878
00879
00880
00881 template <typename TypeTag>
00882 typename MultisegmentWell<TypeTag>::EvalWell
00883 MultisegmentWell<TypeTag>::
00884 volumeFraction(const int seg, const int comp_idx) const
00885 {
00886 const PhaseUsage& pu = phaseUsage();
00887
00888 if (active()[Water] && comp_idx == pu.phase_pos[Water]) {
00889 return primary_variables_evaluation_[seg][WFrac];
00890 }
00891
00892 if (active()[Gas] && comp_idx == pu.phase_pos[Gas]) {
00893 return primary_variables_evaluation_[seg][GFrac];
00894 }
00895
00896
00897 EvalWell oil_fraction = 1.0;
00898 if (active()[Water]) {
00899 oil_fraction -= primary_variables_evaluation_[seg][WFrac];
00900 }
00901
00902 if (active()[Gas]) {
00903 oil_fraction -= primary_variables_evaluation_[seg][GFrac];
00904 }
00905
00906
00907
00908 return oil_fraction;
00909 }
00910
00911
00912
00913
00914
00915 template <typename TypeTag>
00916 typename MultisegmentWell<TypeTag>::EvalWell
00917 MultisegmentWell<TypeTag>::
00918 volumeFractionScaled(const int seg, const int comp_idx) const
00919 {
00920
00921
00922
00923 const double scale = scalingFactor(comp_idx);
00924 if (scale > 0.) {
00925 return volumeFraction(seg, comp_idx) / scale;
00926 }
00927
00928 return volumeFraction(seg, comp_idx);
00929 }
00930
00931
00932
00933
00934
00935 template <typename TypeTag>
00936 typename MultisegmentWell<TypeTag>::EvalWell
00937 MultisegmentWell<TypeTag>::
00938 surfaceVolumeFraction(const int seg, const int comp_idx) const
00939 {
00940 EvalWell sum_volume_fraction_scaled = 0.;
00941 const int num_comp = numComponents();
00942 for (int idx = 0; idx < num_comp; ++idx) {
00943 sum_volume_fraction_scaled += volumeFractionScaled(seg, idx);
00944 }
00945
00946 assert(sum_volume_fraction_scaled.value() != 0.);
00947
00948 return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
00949 }
00950
00951
00952
00953
00954
00955 template <typename TypeTag>
00956 void
00957 MultisegmentWell<TypeTag>::
00958 computePerfRate(const IntensiveQuantities& int_quants,
00959 const std::vector<EvalWell>& mob_perfcells,
00960 const int seg,
00961 const int perf,
00962 const EvalWell& segment_pressure,
00963 const bool& allow_cf,
00964 std::vector<EvalWell>& cq_s) const
00965 {
00966 const int num_comp = numComponents();
00967 std::vector<EvalWell> cmix_s(num_comp, 0.0);
00968
00969
00970 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
00971 cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
00972 }
00973
00974 const auto& fs = int_quants.fluidState();
00975
00976 const EvalWell pressure_cell = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
00977 const EvalWell rs = extendEval(fs.Rs());
00978 const EvalWell rv = extendEval(fs.Rv());
00979
00980
00981 std::vector<EvalWell> b_perfcells(num_comp, 0.0);
00982
00983 for (int phase = 0; phase < number_of_phases_; ++phase) {
00984 const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase);
00985 b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos));
00986 }
00987
00988
00989 const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
00990
00991 const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];
00992
00993
00994
00995 const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff);
00996
00997 const Opm::PhaseUsage& pu = phaseUsage();
00998
00999
01000 if ( drawdown > 0.0) {
01001
01002 if (!allow_cf && well_type_ == INJECTOR) {
01003 return;
01004 }
01005
01006
01007 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01008 const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown);
01009 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
01010 }
01011
01012 if (active()[Oil] && active()[Gas]) {
01013 const int oilpos = pu.phase_pos[Oil];
01014 const int gaspos = pu.phase_pos[Gas];
01015 const EvalWell cq_s_oil = cq_s[oilpos];
01016 const EvalWell cq_s_gas = cq_s[gaspos];
01017 cq_s[gaspos] += rs * cq_s_oil;
01018 cq_s[oilpos] += rv * cq_s_gas;
01019 }
01020 } else {
01021
01022 if (!allow_cf && well_type_ == PRODUCER) {
01023 return;
01024 }
01025
01026
01027 EvalWell total_mob = mob_perfcells[0];
01028 for (int comp_idx = 1; comp_idx < num_comp; ++comp_idx) {
01029 total_mob += mob_perfcells[comp_idx];
01030 }
01031
01032
01033 const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown);
01034
01035
01036 EvalWell volume_ratio = 0.0;
01037 if (active()[Water]) {
01038 const int watpos = pu.phase_pos[Water];
01039 volume_ratio += cmix_s[watpos] / b_perfcells[watpos];
01040 }
01041
01042 if (active()[Oil] && active()[Gas]) {
01043 const int oilpos = pu.phase_pos[Oil];
01044 const int gaspos = pu.phase_pos[Gas];
01045
01046
01047
01048
01049 const EvalWell d = 1.0 - rv * rs;
01050
01051 if (d.value() == 0.0) {
01052 OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << name() << " during flux calcuation"
01053 << " with rs " << rs << " and rv " << rv);
01054 }
01055
01056 const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
01057 volume_ratio += tmp_oil / b_perfcells[oilpos];
01058
01059 const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
01060 volume_ratio += tmp_gas / b_perfcells[gaspos];
01061 } else {
01062 if (active()[Oil]) {
01063 const int oilpos = pu.phase_pos[Oil];
01064 volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos];
01065 }
01066 if (active()[Gas]) {
01067 const int gaspos = pu.phase_pos[Gas];
01068 volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos];
01069 }
01070 }
01071
01072 EvalWell cqt_is = cqt_i / volume_ratio;
01073 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01074 cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
01075 }
01076 }
01077 }
01078
01079
01080
01081
01082
01083 template <typename TypeTag>
01084 typename MultisegmentWell<TypeTag>::EvalWell
01085 MultisegmentWell<TypeTag>::
01086 extendEval(const Eval& in) const
01087 {
01088 EvalWell out = 0.0;
01089 out.setValue(in.value());
01090 for(int eq_idx = 0; eq_idx < numEq;++eq_idx) {
01091 out.setDerivative(eq_idx, in.derivative(eq_idx));
01092 }
01093 return out;
01094 }
01095
01096
01097
01098
01099
01100 template <typename TypeTag>
01101 void
01102 MultisegmentWell<TypeTag>::
01103 computeSegmentFluidProperties(const Simulator& ebosSimulator)
01104 {
01105
01106
01107
01108
01109
01110
01111
01112 EvalWell temperature;
01113
01114
01115
01116
01117
01118 int pvt_region_index;
01119 {
01120
01121 const int cell_idx = well_cells_[0];
01122 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
01123 const auto& fs = intQuants.fluidState();
01124 temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
01125 pvt_region_index = fs.pvtRegionIndex();
01126 }
01127
01128 std::vector<double> surf_dens(number_of_phases_);
01129
01130
01131 for (int phase = 0; phase < number_of_phases_; ++phase) {
01132 surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index );
01133 }
01134
01135 const int num_comp = numComponents();
01136 const Opm::PhaseUsage& pu = phaseUsage();
01137 for (int seg = 0; seg < numberOfSegments(); ++seg) {
01138
01139 std::vector<EvalWell> mix_s(num_comp, 0.0);
01140 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01141 mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
01142 }
01143
01144 std::vector<EvalWell> b(num_comp, 0.0);
01145
01146 std::vector<EvalWell> visc(number_of_phases_, 0.0);
01147 const EvalWell seg_pressure = getSegmentPressure(seg);
01148 if (pu.phase_used[BlackoilPhases::Aqua]) {
01149
01150 const int water_pos = pu.phase_pos[BlackoilPhases::Aqua];
01151 b[water_pos] =
01152 FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
01153 visc[water_pos] =
01154 FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure);
01155 }
01156
01157 EvalWell rv(0.0);
01158
01159 if (pu.phase_used[BlackoilPhases::Vapour]) {
01160 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
01161 if (pu.phase_used[BlackoilPhases::Liquid]) {
01162 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
01163 const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
01164 if (mix_s[oilpos] > 0.0) {
01165 if (mix_s[gaspos] > 0.0) {
01166 rv = mix_s[oilpos] / mix_s[gaspos];
01167 }
01168
01169 if (rv > rvmax) {
01170 rv = rvmax;
01171 }
01172 b[gaspos] =
01173 FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
01174 visc[gaspos] =
01175 FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
01176 } else {
01177 b[gaspos] =
01178 FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
01179 visc[gaspos] =
01180 FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
01181 }
01182 } else {
01183
01184 b[gaspos] =
01185 FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
01186 visc[gaspos] =
01187 FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
01188 }
01189 }
01190
01191 EvalWell rs(0.0);
01192
01193 if (pu.phase_used[BlackoilPhases::Liquid]) {
01194 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
01195 if (pu.phase_used[BlackoilPhases::Liquid]) {
01196 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
01197 const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
01198 if (mix_s[gaspos] > 0.0) {
01199 if (mix_s[oilpos] > 0.0) {
01200 rs = mix_s[gaspos] / mix_s[oilpos];
01201 }
01202
01203 if (rs > rsmax) {
01204 rs = rsmax;
01205 }
01206 b[oilpos] =
01207 FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
01208 visc[oilpos] =
01209 FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
01210 } else {
01211 b[oilpos] =
01212 FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
01213 visc[oilpos] =
01214 FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
01215 }
01216 } else {
01217
01218 b[oilpos] =
01219 FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
01220 visc[oilpos] =
01221 FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
01222 }
01223 }
01224
01225 std::vector<EvalWell> mix(mix_s);
01226 if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
01227 const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
01228 const int oilpos = pu.phase_pos[BlackoilPhases::Liquid];
01229 if (rs != 0.0) {
01230 mix[gaspos] = (mix_s[gaspos] - mix_s[oilpos] * rs) / (1. - rs * rv);
01231 }
01232 if (rv != 0.0) {
01233 mix[oilpos] = (mix_s[oilpos] - mix_s[gaspos] * rv) / (1. - rs * rv);
01234 }
01235 }
01236
01237 EvalWell volrat(0.0);
01238 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01239 volrat += mix[comp_idx] / b[comp_idx];
01240 }
01241
01242 segment_viscosities_[seg] = 0.;
01243
01244 for (int p = 0; p < number_of_phases_; ++p) {
01245 const EvalWell phase_fraction = mix[p] / b[p] / volrat;
01246 segment_viscosities_[seg] += visc[p] * phase_fraction;
01247 }
01248
01249 EvalWell density(0.0);
01250 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01251 density += surf_dens[comp_idx] * mix_s[comp_idx];
01252 }
01253 segment_densities_[seg] = density / volrat;
01254
01255
01256 segment_mass_rates_[seg] = 0.;
01257 for (int phase = 0; phase < number_of_phases_; ++phase) {
01258 const EvalWell rate = getSegmentRate(seg, phase);
01259 segment_mass_rates_[seg] += rate * surf_dens[phase];
01260 }
01261 }
01262 }
01263
01264
01265
01266
01267
01268 template <typename TypeTag>
01269 typename MultisegmentWell<TypeTag>::EvalWell
01270 MultisegmentWell<TypeTag>::
01271 getSegmentPressure(const int seg) const
01272 {
01273 return primary_variables_evaluation_[seg][SPres];
01274 }
01275
01276
01277
01278
01279
01280 template <typename TypeTag>
01281 typename MultisegmentWell<TypeTag>::EvalWell
01282 MultisegmentWell<TypeTag>::
01283 getSegmentRate(const int seg,
01284 const int comp_idx) const
01285 {
01286 return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx);
01287 }
01288
01289
01290
01291
01292
01293 template <typename TypeTag>
01294 typename MultisegmentWell<TypeTag>::EvalWell
01295 MultisegmentWell<TypeTag>::
01296 getSegmentGTotal(const int seg) const
01297 {
01298 return primary_variables_evaluation_[seg][GTotal];
01299 }
01300
01301
01302
01303
01304
01305 template <typename TypeTag>
01306 void
01307 MultisegmentWell<TypeTag>::
01308 getMobility(const Simulator& ebosSimulator,
01309 const int perf,
01310 std::vector<EvalWell>& mob) const
01311 {
01312
01313 const int np = number_of_phases_;
01314 const int cell_idx = well_cells_[perf];
01315 assert (int(mob.size()) == numComponents());
01316 const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
01317 const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
01318
01319
01320
01321 const int satid = saturation_table_number_[perf] - 1;
01322 const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
01323 if( satid == satid_elem ) {
01324
01325 for (int phase = 0; phase < np; ++phase) {
01326 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
01327 mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
01328 }
01329
01330
01331
01332 } else {
01333
01334 const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
01335 Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
01336 MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
01337
01338
01339 materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
01340
01341
01342 for (int phase = 0; phase < np; ++phase) {
01343 int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
01344 mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
01345 }
01346
01347 }
01348 }
01349
01350
01351
01352
01353
01354 template <typename TypeTag>
01355 void
01356 MultisegmentWell<TypeTag>::
01357 assembleControlEq() const
01358 {
01359 EvalWell control_eq(0.0);
01360
01361 switch (well_controls_get_current_type(well_controls_)) {
01362 case THP:
01363 {
01364 OPM_THROW(std::runtime_error, "Not handling THP control for Multisegment wells for now");
01365 }
01366 case BHP:
01367 {
01368 const double target_bhp = well_controls_get_current_target(well_controls_);
01369 control_eq = getSegmentPressure(0) - target_bhp;
01370 break;
01371 }
01372 case SURFACE_RATE:
01373 {
01374
01375 int number_phases_under_control = 0;
01376 const double* distr = well_controls_get_current_distr(well_controls_);
01377 for (int phase = 0; phase < number_of_phases_; ++phase) {
01378 if (distr[phase] > 0.0) {
01379 ++number_phases_under_control;
01380 }
01381 }
01382 assert(number_phases_under_control > 0);
01383
01384 const std::vector<double> g = {1.0, 1.0, 0.01};
01385 const double target_rate = well_controls_get_current_target(well_controls_);
01386
01387
01388 if (number_phases_under_control == 1) {
01389 for (int phase = 0; phase < number_of_phases_; ++phase) {
01390 if (distr[phase] > 0.) {
01391 control_eq = getSegmentGTotal(0) * volumeFraction(0, phase) - g[phase] * target_rate;
01392 break;
01393 }
01394 }
01395 } else {
01396 EvalWell rate_for_control(0.0);
01397 const EvalWell G_total = getSegmentGTotal(0);
01398 for (int phase = 0; phase < number_of_phases_; ++phase) {
01399 if (distr[phase] > 0.) {
01400 rate_for_control += G_total * volumeFractionScaled(0, phase);
01401 }
01402 }
01403
01404 control_eq = rate_for_control - target_rate;
01405 }
01406 break;
01407 }
01408 case RESERVOIR_RATE:
01409 {
01410 EvalWell rate_for_control(0.0);
01411 const double* distr = well_controls_get_current_distr(well_controls_);
01412 for (int phase = 0; phase < number_of_phases_; ++phase) {
01413 if (distr[phase] > 0.) {
01414 rate_for_control += getSegmentGTotal(0) * volumeFraction(0, phase);
01415 }
01416 }
01417 const double target_rate = well_controls_get_current_target(well_controls_);
01418 control_eq = rate_for_control - target_rate;
01419 break;
01420 }
01421 default:
01422 OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name());
01423 }
01424
01425
01426
01427
01428 resWell_[0][SPres] = control_eq.value();
01429 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01430 duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq);
01431 }
01432 }
01433
01434
01435
01436
01437
01438 template <typename TypeTag>
01439 void
01440 MultisegmentWell<TypeTag>::
01441 assemblePressureEq(const int seg) const
01442 {
01443 assert(seg != 0);
01444
01445
01446 EvalWell pressure_equation = getSegmentPressure(seg);
01447
01448
01449
01450 pressure_equation -= getHydroPressureLoss(seg);
01451
01452 if (frictionalPressureLossConsidered()) {
01453 pressure_equation -= getFrictionPressureLoss(seg);
01454 }
01455
01456 resWell_[seg][SPres] = pressure_equation.value();
01457 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01458 duneD_[seg][seg][SPres][pv_idx] = pressure_equation.derivative(pv_idx + numEq);
01459 }
01460
01461
01462 const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
01463 const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
01464
01465 resWell_[seg][SPres] -= outlet_pressure.value();
01466 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01467 duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
01468 }
01469
01470 if (accelerationalPressureLossConsidered()) {
01471 handleAccelerationPressureLoss(seg);
01472 }
01473 }
01474
01475
01476
01477
01478
01479 template <typename TypeTag>
01480 typename MultisegmentWell<TypeTag>::EvalWell
01481 MultisegmentWell<TypeTag>::
01482 getHydroPressureLoss(const int seg) const
01483 {
01484 return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg];
01485 }
01486
01487
01488
01489
01490
01491 template <typename TypeTag>
01492 typename MultisegmentWell<TypeTag>::EvalWell
01493 MultisegmentWell<TypeTag>::
01494 getFrictionPressureLoss(const int seg) const
01495 {
01496 const EvalWell mass_rate = segment_mass_rates_[seg];
01497 const EvalWell density = segment_densities_[seg];
01498 const EvalWell visc = segment_viscosities_[seg];
01499 const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
01500 const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength();
01501 assert(length > 0.);
01502 const double roughness = segmentSet()[seg].roughness();
01503 const double area = segmentSet()[seg].crossArea();
01504 const double diameter = segmentSet()[seg].internalDiameter();
01505
01506 const double sign = mass_rate < 0. ? 1.0 : - 1.0;
01507
01508 return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
01509 }
01510
01511
01512
01513
01514
01515 template <typename TypeTag>
01516 void
01517 MultisegmentWell<TypeTag>::
01518 handleAccelerationPressureLoss(const int seg) const
01519 {
01520
01521
01522 const double area = segmentSet()[seg].crossArea();
01523 const EvalWell mass_rate = segment_mass_rates_[seg];
01524 const EvalWell density = segment_densities_[seg];
01525 const EvalWell out_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
01526
01527 resWell_[seg][SPres] -= out_velocity_head.value();
01528 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01529 duneD_[seg][seg][SPres][pv_idx] -= out_velocity_head.derivative(pv_idx + numEq);
01530 }
01531
01532
01533 double max_area = area;
01534 for (const int inlet : segment_inlets_[seg]) {
01535 const double inlet_area = segmentSet()[inlet].crossArea();
01536 if (inlet_area > max_area) {
01537 max_area = inlet_area;
01538 }
01539 }
01540
01541
01542 for (const int inlet : segment_inlets_[seg]) {
01543 const EvalWell density = segment_densities_[inlet];
01544 const EvalWell mass_rate = segment_mass_rates_[inlet];
01545 const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
01546 resWell_[seg][SPres] += inlet_velocity_head.value();
01547 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01548 duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq);
01549 }
01550 }
01551 }
01552
01553
01554
01555
01556
01557 template <typename TypeTag>
01558 void
01559 MultisegmentWell<TypeTag>::
01560 processFractions(const int seg) const
01561 {
01562 const PhaseUsage& pu = phaseUsage();
01563
01564 std::vector<double> fractions(number_of_phases_, 0.0);
01565
01566 assert( active()[Oil] );
01567 const int oil_pos = pu.phase_pos[Oil];
01568 fractions[oil_pos] = 1.0;
01569
01570 if ( active()[Water] ) {
01571 const int water_pos = pu.phase_pos[Water];
01572 fractions[water_pos] = primary_variables_[seg][WFrac];
01573 fractions[oil_pos] -= fractions[water_pos];
01574 }
01575
01576 if ( active()[Gas] ) {
01577 const int gas_pos = pu.phase_pos[Gas];
01578 fractions[gas_pos] = primary_variables_[seg][GFrac];
01579 fractions[oil_pos] -= fractions[gas_pos];
01580 }
01581
01582 if (active()[Water]) {
01583 const int water_pos = pu.phase_pos[Water];
01584 if (fractions[water_pos] < 0.0) {
01585 if ( active()[Gas] ) {
01586 fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
01587 }
01588 fractions[oil_pos] /= (1.0 - fractions[water_pos]);
01589 fractions[water_pos] = 0.0;
01590 }
01591 }
01592
01593 if (active()[Gas]) {
01594 const int gas_pos = pu.phase_pos[Gas];
01595 if (fractions[gas_pos] < 0.0) {
01596 if ( active()[Water] ) {
01597 fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
01598 }
01599 fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
01600 fractions[gas_pos] = 0.0;
01601 }
01602 }
01603
01604 if (fractions[oil_pos] < 0.0) {
01605 if ( active()[Water] ) {
01606 fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
01607 }
01608 if ( active()[Gas] ) {
01609 fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
01610 }
01611 fractions[oil_pos] = 0.0;
01612 }
01613
01614 if ( active()[Water] ) {
01615 primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
01616 }
01617
01618 if ( active()[Gas] ) {
01619 primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
01620 }
01621 }
01622
01623
01624
01625
01626
01627 template <typename TypeTag>
01628 void
01629 MultisegmentWell<TypeTag>::
01630 updateWellStateFromPrimaryVariables(WellState& well_state) const
01631 {
01632 const PhaseUsage& pu = phaseUsage();
01633 assert( active()[Oil] );
01634 const int oil_pos = pu.phase_pos[Oil];
01635
01636 for (int seg = 0; seg < numberOfSegments(); ++seg) {
01637 std::vector<double> fractions(number_of_phases_, 0.0);
01638 fractions[oil_pos] = 1.0;
01639
01640 if ( active()[Water] ) {
01641 const int water_pos = pu.phase_pos[Water];
01642 fractions[water_pos] = primary_variables_[seg][WFrac];
01643 fractions[oil_pos] -= fractions[water_pos];
01644 }
01645
01646 if ( active()[Gas] ) {
01647 const int gas_pos = pu.phase_pos[Gas];
01648 fractions[gas_pos] = primary_variables_[seg][GFrac];
01649 fractions[oil_pos] -= fractions[gas_pos];
01650 }
01651
01652
01653 for (int p = 0; p < number_of_phases_; ++p) {
01654 const double scale = scalingFactor(p);
01655
01656 if (scale > 0.) {
01657 fractions[p] /= scale;
01658 } else {
01659
01660 fractions[p] = 0.;
01661 }
01662 }
01663
01664
01665 const double g_total = primary_variables_[seg][GTotal];
01666 const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
01667 for (int p = 0; p < number_of_phases_; ++p) {
01668 const double phase_rate = g_total * fractions[p];
01669 well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate;
01670 if (seg == 0) {
01671 well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate;
01672 }
01673 }
01674
01675
01676 well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres];
01677 if (seg == 0) {
01678 well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index];
01679 }
01680 }
01681 }
01682
01683
01684
01685
01686
01687 template <typename TypeTag>
01688 double
01689 MultisegmentWell<TypeTag>::
01690 scalingFactor(const int comp_idx) const
01691 {
01692 const double* distr = well_controls_get_current_distr(well_controls_);
01693
01694 if (well_controls_get_current_type(well_controls_) == RESERVOIR_RATE) {
01695
01696
01697 return distr[comp_idx];
01698 }
01699
01700 const PhaseUsage& pu = phaseUsage();
01701
01702 if (active()[Water] && pu.phase_pos[Water] == comp_idx)
01703 return 1.0;
01704 if (active()[Oil] && pu.phase_pos[Oil] == comp_idx)
01705 return 1.0;
01706 if (active()[Gas] && pu.phase_pos[Gas] == comp_idx)
01707 return 0.01;
01708
01709
01710
01711
01712 assert(false);
01713 return 1.0;
01714 }
01715
01716
01717
01718
01719
01720 template <typename TypeTag>
01721 bool
01722 MultisegmentWell<TypeTag>::
01723 frictionalPressureLossConsidered() const
01724 {
01725
01726 return (segmentSet().compPressureDrop() != WellSegment::H__);
01727 }
01728
01729
01730
01731
01732
01733 template <typename TypeTag>
01734 bool
01735 MultisegmentWell<TypeTag>::
01736 accelerationalPressureLossConsidered() const
01737 {
01738 return (segmentSet().compPressureDrop() == WellSegment::HFA);
01739 }
01740
01741
01742
01743
01744
01745 template<typename TypeTag>
01746 void
01747 MultisegmentWell<TypeTag>::
01748 iterateWellEquations(Simulator& ebosSimulator,
01749 const double dt,
01750 WellState& well_state)
01751 {
01752
01753
01754
01755
01756
01757 const int max_iter_number = param_.max_inner_iter_ms_wells_;
01758 int it = 0;
01759 for (; it < max_iter_number; ++it) {
01760
01761 assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true);
01762
01763 #if HAVE_UMFPACK
01764 const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
01765 #else
01766 const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
01767 #endif
01768
01769
01770
01771
01772
01773
01774 const std::vector<double> B {0.5, 0.5, 0.005};
01775
01776 const ConvergenceReport report = getWellConvergence(B);
01777
01778 if (report.converged) {
01779 break;
01780 }
01781
01782 updateWellState(dx_well, true, well_state);
01783
01784 initPrimaryVariablesEvaluation();
01785 }
01786
01787 }
01788
01789
01790
01791
01792
01793 template<typename TypeTag>
01794 void
01795 MultisegmentWell<TypeTag>::
01796 assembleWellEqWithoutIteration(Simulator& ebosSimulator,
01797 const double dt,
01798 WellState& well_state,
01799 bool only_wells)
01800 {
01801
01802 computeSegmentFluidProperties(ebosSimulator);
01803
01804
01805 if (!only_wells) {
01806 duneB_ = 0.0;
01807 duneC_ = 0.0;
01808 }
01809
01810 duneD_ = 0.0;
01811 resWell_ = 0.0;
01812
01813
01814
01815
01816
01817
01818 auto& ebosJac = ebosSimulator.model().linearizer().matrix();
01819 auto& ebosResid = ebosSimulator.model().linearizer().residual();
01820
01821 const bool allow_cf = getAllowCrossFlow();
01822
01823 const int nseg = numberOfSegments();
01824 const int num_comp = numComponents();
01825
01826 for (int seg = 0; seg < nseg; ++seg) {
01827
01828
01829 {
01830 const double volume = segmentSet()[seg].volume();
01831
01832 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01833 EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx])
01834 + getSegmentRate(seg, comp_idx);
01835
01836 resWell_[seg][comp_idx] += accumulation_term.value();
01837 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01838 duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
01839 }
01840 }
01841 }
01842
01843
01844 {
01845 for (const int inlet : segment_inlets_[seg]) {
01846 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01847 const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx);
01848 resWell_[seg][comp_idx] -= inlet_rate.value();
01849 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01850 duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq);
01851 }
01852 }
01853 }
01854 }
01855
01856
01857 const EvalWell seg_pressure = getSegmentPressure(seg);
01858 for (const int perf : segment_perforations_[seg]) {
01859 const int cell_idx = well_cells_[perf];
01860 const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0));
01861 std::vector<EvalWell> mob(num_comp, 0.0);
01862 getMobility(ebosSimulator, perf, mob);
01863 std::vector<EvalWell> cq_s(num_comp, 0.0);
01864 computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s);
01865
01866 for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
01867
01868 const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
01869
01870 if (!only_wells) {
01871
01872
01873
01874
01875 ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
01876 }
01877
01878
01879 resWell_[seg][comp_idx] -= cq_s_effective.value();
01880
01881
01882 for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
01883 if (!only_wells) {
01884
01885 duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq);
01886 }
01887
01888 duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
01889 }
01890
01891 for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
01892 if (!only_wells) {
01893
01894 ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx);
01895 duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
01896 }
01897 }
01898 }
01899
01900
01901
01902 }
01903
01904
01905 if (seg == 0) {
01906 assembleControlEq();
01907 } else {
01908 assemblePressureEq(seg);
01909 }
01910 }
01911 }
01912
01913 }