00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
00025 #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
00026
00027 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
00028
00029 #include <opm/autodiff/AutoDiffBlock.hpp>
00030 #include <opm/autodiff/AutoDiffHelpers.hpp>
00031 #include <opm/autodiff/GridHelpers.hpp>
00032 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00033 #include <opm/autodiff/GeoProps.hpp>
00034 #include <opm/autodiff/WellDensitySegmented.hpp>
00035
00036 #include <opm/core/grid.h>
00037 #include <opm/core/linalg/LinearSolverInterface.hpp>
00038 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00039 #include <opm/core/props/rock/RockCompressibility.hpp>
00040 #include <opm/common/ErrorMacros.hpp>
00041 #include <opm/common/Exceptions.hpp>
00042 #include <opm/parser/eclipse/Units/Units.hpp>
00043 #include <opm/core/well_controls.h>
00044 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00045
00046 #include <cassert>
00047 #include <cmath>
00048 #include <iostream>
00049 #include <iomanip>
00050 #include <limits>
00051
00052 namespace Opm {
00053
00054
00055
00056 namespace detail {
00057
00058 template <class PU>
00059 int polymerPos(const PU& pu)
00060 {
00061 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
00062 int pos = 0;
00063 for (int phase = 0; phase < maxnp; ++phase) {
00064 if (pu.phase_used[phase]) {
00065 pos++;
00066 }
00067 }
00068
00069 return pos;
00070 }
00071
00072 }
00073
00074
00075
00076 template <class Grid>
00077 BlackoilPolymerModel<Grid>::BlackoilPolymerModel(const typename Base::ModelParameters& param,
00078 const Grid& grid,
00079 const BlackoilPropsAdFromDeck& fluid,
00080 const DerivedGeology& geo,
00081 const RockCompressibility* rock_comp_props,
00082 const PolymerPropsAd& polymer_props_ad,
00083 const StandardWells& well_model,
00084 const NewtonIterationBlackoilInterface& linsolver,
00085 std::shared_ptr< const EclipseState > eclipse_state,
00086 const bool has_disgas,
00087 const bool has_vapoil,
00088 const bool has_polymer,
00089 const bool has_plyshlog,
00090 const bool has_shrate,
00091 const std::vector<double>& wells_rep_radius,
00092 const std::vector<double>& wells_perf_length,
00093 const std::vector<double>& wells_bore_diameter,
00094 const bool terminal_output)
00095 : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver, eclipse_state,
00096 has_disgas, has_vapoil, terminal_output),
00097 polymer_props_ad_(polymer_props_ad),
00098 has_polymer_(has_polymer),
00099 has_plyshlog_(has_plyshlog),
00100 has_shrate_(has_shrate),
00101 poly_pos_(detail::polymerPos(fluid.phaseUsage())),
00102 wells_rep_radius_(wells_rep_radius),
00103 wells_perf_length_(wells_perf_length),
00104 wells_bore_diameter_(wells_bore_diameter)
00105 {
00106 if (has_polymer_) {
00107 if (!active_[Water]) {
00108 OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
00109 }
00110 residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169);
00111
00112 sd_.rq.resize(fluid_.numPhases() + 1);
00113 residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
00114 Base::material_name_.push_back("Polymer");
00115 assert(poly_pos_ == fluid_.numPhases());
00116 }
00117 }
00118
00119
00120
00121 template <class Grid>
00122 void
00123 BlackoilPolymerModel<Grid>::
00124 prepareStep(const SimulatorTimerInterface& timer,
00125 const ReservoirState& reservoir_state,
00126 const WellState& well_state)
00127 {
00128 Base::prepareStep(timer, reservoir_state, well_state);
00129 auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
00130
00131
00132 cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
00133 }
00134
00135
00136
00137
00138 template <class Grid>
00139 void
00140 BlackoilPolymerModel<Grid>::
00141 afterStep(const SimulatorTimerInterface& ,
00142 ReservoirState& reservoir_state,
00143 WellState& )
00144 {
00145 computeCmax(reservoir_state);
00146 }
00147
00148
00149
00150
00151
00152 template <class Grid>
00153 void
00154 BlackoilPolymerModel<Grid>::makeConstantState(SolutionState& state) const
00155 {
00156 Base::makeConstantState(state);
00157 state.concentration = ADB::constant(state.concentration.value());
00158 }
00159
00160
00161
00162
00163
00164 template <class Grid>
00165 std::vector<V>
00166 BlackoilPolymerModel<Grid>::variableStateInitials(const ReservoirState& x,
00167 const WellState& xw) const
00168 {
00169 std::vector<V> vars0 = Base::variableStateInitials(x, xw);
00170 assert(int(vars0.size()) == fluid_.numPhases() + 2);
00171
00172
00173 if (has_polymer_) {
00174 const auto& concentration = x.getCellData( x.CONCENTRATION );
00175 assert (not concentration.empty());
00176 const int nc = concentration.size();
00177 const V c = Eigen::Map<const V>(concentration.data() , nc);
00178
00179 auto concentration_pos = vars0.begin() + fluid_.numPhases();
00180 assert(concentration_pos == vars0.end() - 2);
00181 vars0.insert(concentration_pos, c);
00182 }
00183 return vars0;
00184 }
00185
00186
00187
00188
00189
00190 template <class Grid>
00191 std::vector<int>
00192 BlackoilPolymerModel<Grid>::variableStateIndices() const
00193 {
00194 std::vector<int> ind = Base::variableStateIndices();
00195 assert(ind.size() == 5);
00196 if (has_polymer_) {
00197 ind.resize(6);
00198
00199 ind[Concentration] = fluid_.numPhases();
00200
00201 ++ind[Qs];
00202 ++ind[Bhp];
00203 }
00204 return ind;
00205 }
00206
00207
00208
00209
00210 template <class Grid>
00211 typename BlackoilPolymerModel<Grid>::SolutionState
00212 BlackoilPolymerModel<Grid>::variableStateExtractVars(const ReservoirState& x,
00213 const std::vector<int>& indices,
00214 std::vector<ADB>& vars) const
00215 {
00216 SolutionState state = Base::variableStateExtractVars(x, indices, vars);
00217 if (has_polymer_) {
00218 state.concentration = std::move(vars[indices[Concentration]]);
00219 }
00220 return state;
00221 }
00222
00223
00224
00225
00226
00227 template <class Grid>
00228 void
00229 BlackoilPolymerModel<Grid>::computeAccum(const SolutionState& state,
00230 const int aix )
00231 {
00232 Base::computeAccum(state, aix);
00233
00234
00235 if (has_polymer_) {
00236 const ADB& press = state.pressure;
00237 const std::vector<ADB>& sat = state.saturation;
00238 const ADB& c = state.concentration;
00239 const ADB pv_mult = poroMult(press);
00240 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00241
00242 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
00243 const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
00244 const double rho_rock = polymer_props_ad_.rockDensity();
00245 const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
00246 const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
00247
00248 sd_.rq[poly_pos_].accum[aix] = pv_mult * sd_.rq[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
00249 + pv_mult * rho_rock * (1. - phi) / phi * ads;
00250 }
00251
00252 }
00253
00254
00255
00256
00257
00258
00259 template <class Grid>
00260 void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
00261 {
00262 auto& max_concentration = state.getCellData( state.CMAX );
00263 const auto& concentration = state.getCellData( state.CONCENTRATION );
00264 std::transform( max_concentration.begin() ,
00265 max_concentration.end() ,
00266 concentration.begin() ,
00267 max_concentration.begin() ,
00268 [](double c_max , double c) { return std::max( c_max , c ); });
00269 }
00270
00271
00272
00273
00274
00275 template <class Grid>
00276 void
00277 BlackoilPolymerModel<Grid>::
00278 assembleMassBalanceEq(const SolutionState& state)
00279 {
00280
00281
00282
00283
00284
00285
00286
00287
00288 computeAccum(state, 1);
00289
00290
00291
00292 const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
00293 {
00294 const std::vector<ADB> kr = computeRelPerm(state);
00295 for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
00296 sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
00297 }
00298 }
00299
00300
00301 if (has_plyshlog_) {
00302 std::vector<double> water_vel;
00303 std::vector<double> visc_mult;
00304
00305 computeWaterShearVelocityFaces(transi, state.canonical_phase_pressures, state, water_vel, visc_mult);
00306 if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
00307
00308 OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
00309 }
00310 }
00311
00312 for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
00313 const std::vector<PhasePresence>& cond = phaseCondition();
00314 sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
00315 sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
00316 computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
00317
00318 residual_.material_balance_eq[ phaseIdx ] =
00319 pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
00320 + ops_.div*sd_.rq[phaseIdx].mflux;
00321 }
00322
00323
00324
00325
00326
00327
00328 if (active_[ Oil ] && active_[ Gas ]) {
00329 const int po = fluid_.phaseUsage().phase_pos[ Oil ];
00330 const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
00331
00332 const UpwindSelector<double> upwindOil(grid_, ops_,
00333 sd_.rq[po].dh.value());
00334 const ADB rs_face = upwindOil.select(state.rs);
00335
00336 const UpwindSelector<double> upwindGas(grid_, ops_,
00337 sd_.rq[pg].dh.value());
00338 const ADB rv_face = upwindGas.select(state.rv);
00339
00340 residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
00341 residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
00342
00343
00344
00345 }
00346
00347
00348 if (has_polymer_) {
00349 residual_.material_balance_eq[poly_pos_] = pvdt_ * (sd_.rq[poly_pos_].accum[1] - sd_.rq[poly_pos_].accum[0])
00350 + ops_.div*sd_.rq[poly_pos_].mflux;
00351 }
00352
00353
00354 if (param_.update_equations_scaling_) {
00355 updateEquationsScaling();
00356 }
00357
00358 }
00359
00360
00361
00362
00363
00364 template <class Grid>
00365 void BlackoilPolymerModel<Grid>::updateEquationsScaling()
00366 {
00367 Base::updateEquationsScaling();
00368 if (has_polymer_) {
00369 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
00370 residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
00371 }
00372 }
00373
00374
00375
00376
00377
00378 template <class Grid>
00379 void BlackoilPolymerModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
00380 const SolutionState& state,
00381 WellState& xw)
00382
00383 {
00384 Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
00385
00386
00387 if (has_polymer_) {
00388 const ADB mc = computeMc(state);
00389 const int nc = xw.polymerInflow().size();
00390 const V polyin = Eigen::Map<const V>(xw.polymerInflow().data(), nc);
00391 const int nperf = wells().well_connpos[wells().number_of_wells];
00392 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
00393 const V poly_in_perf = subset(polyin, well_cells);
00394 const V poly_mc_perf = subset(mc.value(), well_cells);
00395 const ADB& cq_s_water = cq_s[fluid_.phaseUsage().phase_pos[Water]];
00396 Selector<double> injector_selector(cq_s_water.value());
00397 const V poly_perf = injector_selector.select(poly_in_perf, poly_mc_perf);
00398 const ADB cq_s_poly = cq_s_water * poly_perf;
00399 residual_.material_balance_eq[poly_pos_] -= superset(cq_s_poly, well_cells, nc);
00400 }
00401 }
00402
00403
00404
00405
00406
00407
00408 template <class Grid>
00409 void BlackoilPolymerModel<Grid>::updateState(const V& dx,
00410 ReservoirState& reservoir_state,
00411 WellState& well_state)
00412 {
00413 if (has_polymer_) {
00414
00415 const int np = fluid_.numPhases();
00416 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00417 const V zero = V::Zero(nc);
00418 const int concentration_start = nc * np;
00419 const V dc = subset(dx, Span(nc, 1, concentration_start));
00420
00421
00422 V modified_dx = V::Zero(dx.size() - nc);
00423 modified_dx.head(concentration_start) = dx.head(concentration_start);
00424 const int tail_len = dx.size() - concentration_start - nc;
00425 modified_dx.tail(tail_len) = dx.tail(tail_len);
00426
00427
00428 Base::updateState(modified_dx, reservoir_state, well_state);
00429
00430 {
00431 auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
00432
00433 const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
00434 const V c = (c_old - dc).max(zero);
00435 std::copy(&c[0], &c[0] + nc, concentration.begin());
00436 }
00437 } else {
00438
00439 Base::updateState(dx, reservoir_state, well_state);
00440 }
00441 }
00442
00443
00444
00445
00446
00447 template <class Grid>
00448 void
00449 BlackoilPolymerModel<Grid>::computeMassFlux(const int actph ,
00450 const V& transi,
00451 const ADB& kr ,
00452 const ADB& mu ,
00453 const ADB& rho ,
00454 const ADB& phasePressure,
00455 const SolutionState& state)
00456 {
00457 Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state);
00458
00459
00460 const int canonicalPhaseIdx = canph_[ actph ];
00461 if (canonicalPhaseIdx == Water) {
00462 if (has_polymer_) {
00463 const ADB tr_mult = transMult(state.pressure);
00464 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
00465 const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
00466 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
00467
00468 sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
00469
00470 const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
00471 sd_.rq[poly_pos_].mob = tr_mult * krw_eff * state.concentration * inv_poly_eff_visc;
00472 sd_.rq[poly_pos_].b = sd_.rq[actph].b;
00473 sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
00474 UpwindSelector<double> upwind(grid_, ops_, sd_.rq[poly_pos_].dh.value());
00475
00476 sd_.rq[poly_pos_].mflux = upwind.select(sd_.rq[poly_pos_].b * sd_.rq[poly_pos_].mob) * (transi * sd_.rq[poly_pos_].dh);
00477
00478 sd_.rq[ actph ].mflux = upwind.select(sd_.rq[actph].b * sd_.rq[actph].mob) * (transi * sd_.rq[actph].dh);
00479
00480
00481 if (has_plyshlog_) {
00482 V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
00483 ADB shear_mult_faces_adb = ADB::constant(shear_mult_faces_v);
00484 sd_.rq[poly_pos_].mflux = sd_.rq[poly_pos_].mflux / shear_mult_faces_adb;
00485 sd_.rq[actph].mflux = sd_.rq[actph].mflux / shear_mult_faces_adb;
00486 }
00487 }
00488 }
00489 }
00490
00491
00492
00493
00494 template <class Grid>
00495 SimulatorReport
00496 BlackoilPolymerModel<Grid>::assemble(const ReservoirState& reservoir_state,
00497 WellState& well_state,
00498 const bool initial_assembly)
00499 {
00500 using namespace Opm::AutoDiffGrid;
00501
00502 SimulatorReport report;
00503
00504
00505
00506
00507 wellModel().updateWellControls(well_state);
00508
00509
00510 SolutionState state = variableState(reservoir_state, well_state);
00511
00512 if (initial_assembly) {
00513
00514 SolutionState state0 = state;
00515 makeConstantState(state0);
00516
00517
00518 computeAccum(state0, 0);
00519
00520 wellModel().computeWellConnectionPressures(state0, well_state);
00521 }
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533 assembleMassBalanceEq(state);
00534
00535 if ( ! wellsActive() ) {
00536 return report;
00537 }
00538
00539 std::vector<ADB> mob_perfcells;
00540 std::vector<ADB> b_perfcells;
00541 wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
00542
00543
00544
00545 computeInjectionMobility(state, mob_perfcells);
00546
00547 if (param_.solve_welleq_initially_ && initial_assembly) {
00548
00549 Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
00550 }
00551
00552 V aliveWells;
00553 std::vector<ADB> cq_s;
00554
00555 wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
00556
00557 if (has_plyshlog_) {
00558 std::vector<double> water_vel_wells;
00559 std::vector<double> visc_mult_wells;
00560
00561 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
00562 computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells);
00563
00564 if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) {
00565 OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells ");
00566 }
00567
00568
00569 V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
00570 ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v);
00571 mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
00572 wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
00573 }
00574
00575 wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
00576 wellModel().addWellFluxEq(cq_s, state, residual_);
00577 addWellContributionToMassBalanceEq(cq_s, state, well_state);
00578 wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
00579
00580 report.converged = true;
00581 return report;
00582 }
00583
00584
00585
00586
00587 template <class Grid>
00588 ADB
00589 BlackoilPolymerModel<Grid>::computeMc(const SolutionState& state) const
00590 {
00591 return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
00592 }
00593
00594
00595
00596
00597 template<class Grid>
00598 void
00599 BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi,
00600 const std::vector<ADB>& phasePressure, const SolutionState& state,
00601 std::vector<double>& water_vel, std::vector<double>& visc_mult)
00602 {
00603
00604 const int phase = fluid_.phaseUsage().phase_pos[Water];
00605
00606 const int canonicalPhaseIdx = canph_[phase];
00607
00608 const std::vector<PhasePresence> cond = phaseCondition();
00609
00610 const ADB tr_mult = transMult(state.pressure);
00611 const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
00612 sd_.rq[phase].mob = tr_mult * sd_.rq[phase].kr / mu;
00613
00614
00615 const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
00616 const ADB rhoavg = ops_.caver * rho;
00617 sd_.rq[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
00618 if (use_threshold_pressure_) {
00619 applyThresholdPressures(sd_.rq[ phase ].dh);
00620 }
00621
00622 const ADB& b = sd_.rq[ phase ].b;
00623 const ADB& mob = sd_.rq[ phase ].mob;
00624 const ADB& dh = sd_.rq[ phase ].dh;
00625 UpwindSelector<double> upwind(grid_, ops_, dh.value());
00626
00627 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
00628 ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
00629 cmax,
00630 sd_.rq[phase].kr);
00631 ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
00632 sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
00633
00634 const V& polymer_conc = state.concentration.value();
00635 V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
00636 V visc_mult_faces = upwind.select(visc_mult_cells);
00637
00638 size_t nface = visc_mult_faces.size();
00639 visc_mult.resize(nface);
00640 std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin());
00641
00642 sd_.rq[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
00643
00644
00645 const auto& b_faces_adb = upwind.select(b);
00646 std::vector<double> b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size());
00647
00648 const auto& internal_faces = ops_.internal_faces;
00649
00650 std::vector<double> internal_face_areas;
00651 internal_face_areas.resize(internal_faces.size());
00652
00653 for (int i = 0; i < internal_faces.size(); ++i) {
00654 internal_face_areas[i] = grid_.face_areas[internal_faces[i]];
00655 }
00656
00657 const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
00658 const ADB phiavg_adb = ops_.caver * phi;
00659
00660 std::vector<double> phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size());
00661
00662 water_vel.resize(nface);
00663 std::copy(sd_.rq[0].mflux.value().data(), sd_.rq[0].mflux.value().data() + nface, water_vel.begin());
00664
00665 for (size_t i = 0; i < nface; ++i) {
00666 water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
00667 }
00668
00669
00670 if (has_shrate_) {
00671
00672
00673 const Opm::PhaseUsage pu = fluid_.phaseUsage();
00674 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
00675 const ADB& sw_upwind_adb = upwind.select(sw);
00676 std::vector<double> sw_upwind(sw_upwind_adb.value().data(), sw_upwind_adb.value().data() + sw_upwind_adb.size());
00677
00678
00679 std::vector<double> perm;
00680 perm.resize(transi.size());
00681
00682 for (int i = 0; i < transi.size(); ++i) {
00683 perm[i] = transi[i] / internal_faces[i];
00684 }
00685
00686
00687 const ADB& krw_adb = upwind.select(krw_eff);
00688 std::vector<double> krw_upwind(krw_adb.value().data(), krw_adb.value().data() + krw_adb.size());
00689
00690 const double& shrate_const = polymer_props_ad_.shrate();
00691
00692 const double epsilon = std::numeric_limits<double>::epsilon();
00693
00694
00695
00696 for (size_t i = 0; i < water_vel.size(); ++i) {
00697
00698
00699 if (std::abs(water_vel[i]) < epsilon) {
00700 continue;
00701 }
00702
00703 water_vel[i] *= shrate_const * std::sqrt(phiavg[i] / (perm[i] * sw_upwind[i] * krw_upwind[i]));
00704
00705 }
00706 }
00707
00708 }
00709
00710
00711
00712
00713 template<class Grid>
00714 void
00715 BlackoilPolymerModel<Grid>::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
00716 std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells)
00717 {
00718 if( ! wellsActive() ) return ;
00719
00720 const int nw = wells().number_of_wells;
00721 const int nperf = wells().well_connpos[nw];
00722 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
00723
00724 water_vel_wells.resize(cq_sw.size());
00725 std::copy(cq_sw.value().data(), cq_sw.value().data() + cq_sw.size(), water_vel_wells.begin());
00726
00727 const V& polymer_conc = state.concentration.value();
00728
00729 V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
00730 V visc_mult_wells_v = subset(visc_mult_cells, well_cells);
00731
00732 visc_mult_wells.resize(visc_mult_wells_v.size());
00733 std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin());
00734
00735 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
00736 ADB b_perfcells = subset(sd_.rq[water_pos].b, well_cells);
00737
00738 const ADB& p_perfcells = subset(state.pressure, well_cells);
00739 const V& cdp = wellModel().wellPerforationPressureDiffs();
00740 const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
00741
00742 const ADB drawdown = p_perfcells - perfpressure;
00743
00744
00745 V selectInjectingPerforations = V::Zero(nperf);
00746 for (int c = 0; c < nperf; ++c) {
00747 if (drawdown.value()[c] < 0) {
00748 selectInjectingPerforations[c] = 1;
00749 }
00750 }
00751
00752
00753 for (size_t i = 0; i < well_cells.size(); ++i) {
00754 if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) {
00755 visc_mult_wells[i] = 1.;
00756 }
00757 if (selectInjectingPerforations[i] == 1) {
00758 if (xw.polymerInflow()[well_cells[i]] == 0.) {
00759 visc_mult_wells[i] = 1.;
00760 } else {
00761
00762 const double c_perf = state.concentration.value()[well_cells[i]];
00763 visc_mult_wells[i] = polymer_props_ad_.viscMult(c_perf);
00764 }
00765 }
00766 }
00767
00768 const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
00769 const ADB phi_wells_adb = subset(phi, well_cells);
00770
00771 std::vector<double> phi_wells(phi_wells_adb.value().data(), phi_wells_adb.value().data() + phi_wells_adb.size());
00772
00773 std::vector<double> b_wells(b_perfcells.value().data(), b_perfcells.value().data() + b_perfcells.size());
00774
00775 for (size_t i = 0; i < water_vel_wells.size(); ++i) {
00776 water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]);
00777
00778
00779 }
00780
00781
00782 if (has_shrate_) {
00783 const double& shrate_const = polymer_props_ad_.shrate();
00784 for (size_t i = 0; i < water_vel_wells.size(); ++i) {
00785 water_vel_wells[i] = shrate_const * water_vel_wells[i] / wells_bore_diameter_[i];
00786 }
00787 }
00788
00789 return;
00790 }
00791
00792
00793
00794
00795
00796 template<class Grid>
00797 void
00798 BlackoilPolymerModel<Grid>::
00799 computeInjectionMobility(const SolutionState& state,
00800 std::vector<ADB>& mob_perfcells)
00801 {
00802
00803 if (has_polymer_ && wellModel().localWellsActive()) {
00804 const std::vector<int> well_cells = wellModel().wellOps().well_cells;
00805 const int nperf = well_cells.size();
00806
00807
00808 const ADB& p_perfcells = subset(state.pressure, well_cells);
00809 const V& cdp = wellModel().wellPerforationPressureDiffs();
00810 const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
00811
00812 const ADB drawdown = p_perfcells - perfpressure;
00813
00814
00815 const ADB c_perfcells = subset(state.concentration, well_cells);
00816
00817
00818
00819 std::vector<int> polymer_inj_cells;
00820 std::vector<int> other_well_cells;
00821
00822 polymer_inj_cells.reserve(nperf);
00823 other_well_cells.reserve(nperf);
00824
00825 for (int c = 0; c < nperf; ++c) {
00826
00827 if (drawdown.value()[c] < 0.0 && c_perfcells.value()[c] > 0.0) {
00828 polymer_inj_cells.push_back(c);
00829 } else {
00830 other_well_cells.push_back(c);
00831 }
00832 }
00833
00834
00835 if ( !polymer_inj_cells.empty() ) {
00836
00837 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
00838 const ADB mu_perfcells = subset(sd_.rq[water_pos].mu, well_cells);
00839
00840 const ADB c_poly_inj_cells = subset(c_perfcells, polymer_inj_cells);
00841 const ADB mu_poly_inj_cells = subset(mu_perfcells, polymer_inj_cells);
00842
00843 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(c_poly_inj_cells, mu_poly_inj_cells.value());
00844 const ADB fully_mixing_visc = polymer_props_ad_.viscMult(c_poly_inj_cells) * mu_poly_inj_cells;
00845
00846
00847 ADB mob_polymer_inj = subset(mob_perfcells[water_pos], polymer_inj_cells);
00848 const ADB mob_others = subset(mob_perfcells[water_pos], other_well_cells);
00849
00850 mob_polymer_inj = mob_polymer_inj / inv_wat_eff_visc / fully_mixing_visc;
00851 mob_perfcells[water_pos] = superset(mob_polymer_inj, polymer_inj_cells, nperf) +
00852 superset(mob_others, other_well_cells, nperf);
00853 }
00854 }
00855 }
00856
00857 }
00858
00859 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED