00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
00021 #define OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
00022
00023 #include <opm/autodiff/BlackoilSolventModel.hpp>
00024
00025 #include <opm/autodiff/AutoDiffBlock.hpp>
00026 #include <opm/autodiff/AutoDiffHelpers.hpp>
00027 #include <opm/autodiff/GridHelpers.hpp>
00028 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00029 #include <opm/autodiff/GeoProps.hpp>
00030 #include <opm/autodiff/WellDensitySegmented.hpp>
00031
00032 #include <opm/core/grid.h>
00033 #include <opm/core/linalg/LinearSolverInterface.hpp>
00034 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00035 #include <opm/core/props/rock/RockCompressibility.hpp>
00036 #include <opm/common/ErrorMacros.hpp>
00037 #include <opm/common/Exceptions.hpp>
00038 #include <opm/parser/eclipse/Units/Units.hpp>
00039 #include <opm/core/well_controls.h>
00040 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00041
00042 #include <cassert>
00043 #include <cmath>
00044 #include <iostream>
00045 #include <iomanip>
00046 #include <limits>
00047
00048 namespace Opm {
00049
00050
00051
00052 namespace detail {
00053
00054 template <class PU>
00055 int solventPos(const PU& pu)
00056 {
00057 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
00058 int pos = 0;
00059 for (int phase = 0; phase < maxnp; ++phase) {
00060 if (pu.phase_used[phase]) {
00061 pos++;
00062 }
00063 }
00064
00065 return pos;
00066 }
00067
00068 }
00069
00070
00071
00072 template <class Grid>
00073 BlackoilSolventModel<Grid>::BlackoilSolventModel(const typename Base::ModelParameters& param,
00074 const Grid& grid,
00075 const BlackoilPropsAdFromDeck& fluid,
00076 const DerivedGeology& geo,
00077 const RockCompressibility* rock_comp_props,
00078 const SolventPropsAdFromDeck& solvent_props,
00079 const StandardWellsSolvent& well_model,
00080 const NewtonIterationBlackoilInterface& linsolver,
00081 std::shared_ptr< const EclipseState > eclState,
00082 const bool has_disgas,
00083 const bool has_vapoil,
00084 const bool terminal_output,
00085 const bool has_solvent,
00086 const bool is_miscible)
00087 : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
00088 eclState, has_disgas, has_vapoil, terminal_output),
00089 has_solvent_(has_solvent),
00090 solvent_pos_(detail::solventPos(fluid.phaseUsage())),
00091 solvent_props_(solvent_props),
00092 is_miscible_(is_miscible)
00093
00094 {
00095 if (has_solvent_) {
00096
00097
00098 sd_.rq.resize(fluid_.numPhases() + 1);
00099 residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
00100 Base::material_name_.push_back("Solvent");
00101 assert(solvent_pos_ == fluid_.numPhases());
00102 residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031);
00103
00104 wellModel().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
00105 }
00106 if (is_miscible_) {
00107 mu_eff_.resize(fluid_.numPhases() + 1, ADB::null());
00108 b_eff_.resize(fluid_.numPhases() + 1, ADB::null());
00109 }
00110 }
00111
00112
00113
00114
00115
00116 template <class Grid>
00117 void
00118 BlackoilSolventModel<Grid>::makeConstantState(SolutionState& state) const
00119 {
00120 Base::makeConstantState(state);
00121 state.solvent_saturation = ADB::constant(state.solvent_saturation.value());
00122 }
00123
00124
00125
00126
00127
00128 template <class Grid>
00129 std::vector<V>
00130 BlackoilSolventModel<Grid>::variableStateInitials(const ReservoirState& x,
00131 const WellState& xw) const
00132 {
00133 std::vector<V> vars0 = Base::variableStateInitials(x, xw);
00134 assert(int(vars0.size()) == fluid_.numPhases() + 2);
00135
00136
00137 if (has_solvent_) {
00138 const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL );
00139 const int nc = solvent_saturation.size();
00140 const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
00141
00142
00143
00144 auto solvent_pos = vars0.begin() + fluid_.numPhases();
00145 assert (not solvent_saturation.empty());
00146 assert(solvent_pos == vars0.end() - 2);
00147 vars0.insert(solvent_pos, ss);
00148 }
00149 return vars0;
00150 }
00151
00152
00153
00154
00155
00156 template <class Grid>
00157 std::vector<int>
00158 BlackoilSolventModel<Grid>::variableStateIndices() const
00159 {
00160 std::vector<int> ind = Base::variableStateIndices();
00161 assert(ind.size() == 5);
00162 if (has_solvent_) {
00163 ind.resize(6);
00164
00165 ind[Solvent] = fluid_.numPhases();
00166
00167 ++ind[Qs];
00168 ++ind[Bhp];
00169 }
00170 return ind;
00171 }
00172
00173
00174
00175
00176 template <class Grid>
00177 typename BlackoilSolventModel<Grid>::SolutionState
00178 BlackoilSolventModel<Grid>::variableStateExtractVars(const ReservoirState& x,
00179 const std::vector<int>& indices,
00180 std::vector<ADB>& vars) const
00181 {
00182
00183
00184
00185
00186 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00187 const Opm::PhaseUsage pu = fluid_.phaseUsage();
00188
00189 SolutionState state(fluid_.numPhases());
00190
00191
00192 state.pressure = std::move(vars[indices[Pressure]]);
00193
00194
00195 const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
00196 state.temperature = ADB::constant(temp);
00197
00198
00199 {
00200 ADB so = ADB::constant(V::Ones(nc, 1));
00201
00202 if (active_[ Water ]) {
00203 state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
00204 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
00205 so -= sw;
00206 }
00207 if (has_solvent_) {
00208 state.solvent_saturation = std::move(vars[indices[Solvent]]);
00209 so -= state.solvent_saturation;
00210 }
00211
00212 if (active_[ Gas ]) {
00213
00214
00215 const ADB& xvar = vars[indices[Xvar]];
00216 ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
00217 sg = Base::isSg_*xvar + Base::isRv_*so;
00218 so -= sg;
00219
00220 if (active_[ Oil ]) {
00221
00222 state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation);
00223 sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
00224 if (has_disgas_) {
00225 state.rs = (1-Base::isRs_)*sd_.rsSat + Base::isRs_*xvar;
00226 } else {
00227 state.rs = sd_.rsSat;
00228 }
00229 sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
00230 if (has_vapoil_) {
00231 state.rv = (1-Base::isRv_)*sd_.rvSat + Base::isRv_*xvar;
00232 } else {
00233 state.rv = sd_.rvSat;
00234 }
00235 }
00236 }
00237
00238 if (active_[ Oil ]) {
00239
00240 state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
00241 }
00242 }
00243
00244 wellModel().variableStateExtractWellsVars(indices, vars, state);
00245 return state;
00246
00247 }
00248
00249
00250
00251
00252
00253 template <class Grid>
00254 void
00255 BlackoilSolventModel<Grid>::computeAccum(const SolutionState& state,
00256 const int aix )
00257 {
00258
00259 if (is_miscible_) {
00260 computeEffectiveProperties(state);
00261 }
00262 Base::computeAccum(state, aix);
00263
00264
00265 if (has_solvent_) {
00266 const ADB& press = state.pressure;
00267 const ADB& ss = state.solvent_saturation;
00268 const ADB pv_mult = poroMult(press);
00269 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00270
00271 const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
00272 const std::vector<PhasePresence>& cond = phaseCondition();
00273 sd_.rq[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
00274 sd_.rq[solvent_pos_].accum[aix] = pv_mult * sd_.rq[solvent_pos_].b * ss;
00275 }
00276 }
00277
00278
00279
00280
00281
00282 template <class Grid>
00283 void
00284 BlackoilSolventModel<Grid>::
00285 assembleMassBalanceEq(const SolutionState& state)
00286 {
00287
00288 Base::assembleMassBalanceEq(state);
00289
00290 if (has_solvent_) {
00291 residual_.material_balance_eq[ solvent_pos_ ] =
00292 pvdt_ * (sd_.rq[solvent_pos_].accum[1] - sd_.rq[solvent_pos_].accum[0])
00293 + ops_.div*sd_.rq[solvent_pos_].mflux;
00294 }
00295
00296 }
00297
00298 template <class Grid>
00299 void
00300 BlackoilSolventModel<Grid>::updateEquationsScaling()
00301 {
00302 Base::updateEquationsScaling();
00303 assert(MaxNumPhases + 1 == residual_.matbalscale.size());
00304 if (has_solvent_) {
00305 const ADB& temp_b = sd_.rq[solvent_pos_].b;
00306 ADB::V B = 1. / temp_b.value();
00307 #if HAVE_MPI
00308 if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00309 {
00310 const ParallelISTLInformation& real_info =
00311 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
00312 double B_global_sum = 0;
00313 real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
00314 residual_.matbalscale[solvent_pos_] = B_global_sum / Base::global_nc_;
00315 }
00316 else
00317 #endif
00318 {
00319 residual_.matbalscale[solvent_pos_] = B.mean();
00320 }
00321 }
00322 }
00323
00324 template <class Grid>
00325 void BlackoilSolventModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
00326 const SolutionState& state,
00327 WellState& xw)
00328
00329 {
00330
00331
00332
00333 Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
00334
00335 if (has_solvent_) {
00336 const int nperf = wells().well_connpos[wells().number_of_wells];
00337 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00338
00339 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00340 const ADB zero = ADB::constant(V::Zero(nc));
00341 const ADB& ss = state.solvent_saturation;
00342 const ADB& sg = (active_[ Gas ]
00343 ? state.saturation[ pu.phase_pos[ Gas ] ]
00344 : zero);
00345
00346 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
00347 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00348 ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
00349
00350 const int nw = wells().number_of_wells;
00351 V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
00352
00353 V isProducer = V::Zero(nperf);
00354 V ones = V::Constant(nperf,1.0);
00355 for (int w = 0; w < nw; ++w) {
00356 if(wells().type[w] == PRODUCER) {
00357 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
00358 isProducer[perf] = 1;
00359 }
00360 }
00361 }
00362
00363 const ADB& rs_perfcells = subset(state.rs, well_cells);
00364 const ADB& rv_perfcells = subset(state.rv, well_cells);
00365 int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
00366 int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
00367
00368
00369 const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]) / (ones - rs_perfcells * rv_perfcells);
00370
00371
00372
00373 residual_.material_balance_eq[solvent_pos_] -= superset(cq_s_solvent, well_cells, nc);
00374
00375
00376
00377 residual_.material_balance_eq[gas_pos] += superset(cq_s_solvent, well_cells, nc);
00378
00379 }
00380 }
00381
00382
00383
00384
00385
00386 template <class Grid>
00387 void
00388 BlackoilSolventModel<Grid>::
00389 updateState(const V& dx,
00390 ReservoirState& reservoir_state,
00391 WellState& well_state)
00392 {
00393
00394
00395
00396
00397
00398 using namespace Opm::AutoDiffGrid;
00399 const int np = fluid_.numPhases();
00400 const int nc = numCells(grid_);
00401 const V null;
00402 assert(null.size() == 0);
00403 const V zero = V::Zero(nc);
00404 const V ones = V::Constant(nc,1.0);
00405
00406
00407 const V dp = subset(dx, Span(nc));
00408 int varstart = nc;
00409 const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
00410 varstart += dsw.size();
00411
00412 const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
00413 varstart += dxvar.size();
00414
00415 const V dss = has_solvent_ ? subset(dx, Span(nc, 1, varstart)) : null;
00416 varstart += dss.size();
00417
00418
00419 const V dwells = subset(dx, Span(wellModel().numWellVars(), 1, varstart));
00420 varstart += dwells.size();
00421
00422 assert(varstart == dx.size());
00423
00424
00425 const double dpmaxrel = dpMaxRel();
00426 const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
00427 const V absdpmax = dpmaxrel*p_old.abs();
00428 const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
00429 const V p = (p_old - dp_limited).max(zero);
00430 std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
00431
00432
00433
00434 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00435 const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
00436 const double dsmax = dsMax();
00437
00438
00439
00440 V so;
00441 V sw = zero;
00442 V sg = zero;
00443 V ss = zero;
00444
00445
00446
00447 {
00448 V maxVal = zero;
00449 V dso = zero;
00450 if (active_[Water]){
00451 maxVal = dsw.abs().max(maxVal);
00452 dso = dso - dsw;
00453 }
00454
00455 V dsg;
00456 if (active_[Gas]){
00457 dsg = Base::isSg_ * dxvar - Base::isRv_ * dsw;
00458 maxVal = dsg.abs().max(maxVal);
00459 dso = dso - dsg;
00460 }
00461 if (has_solvent_){
00462 maxVal = dss.abs().max(maxVal);
00463
00464
00465
00466 }
00467
00468 maxVal = dso.abs().max(maxVal);
00469
00470 V step = dsmax/maxVal;
00471 step = step.min(1.);
00472
00473 if (active_[Water]) {
00474 const int pos = pu.phase_pos[ Water ];
00475 const V sw_old = s_old.col(pos);
00476 sw = sw_old - step * dsw;
00477 }
00478
00479 if (active_[Gas]) {
00480 const int pos = pu.phase_pos[ Gas ];
00481 const V sg_old = s_old.col(pos);
00482 sg = sg_old - step * dsg;
00483 }
00484 if (has_solvent_) {
00485 auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
00486 const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
00487 ss = ss_old - step * dss;
00488 }
00489
00490 const int pos = pu.phase_pos[ Oil ];
00491 const V so_old = s_old.col(pos);
00492 so = so_old - step * dso;
00493 }
00494
00495
00496 auto ixg = sg < 0;
00497 for (int c = 0; c < nc; ++c) {
00498 if (ixg[c]) {
00499 sw[c] = sw[c] / (1-sg[c]);
00500 so[c] = so[c] / (1-sg[c]);
00501 sg[c] = 0;
00502 }
00503 }
00504
00505 auto ixo = so < 0;
00506 for (int c = 0; c < nc; ++c) {
00507 if (ixo[c]) {
00508 sw[c] = sw[c] / (1-so[c]);
00509 sg[c] = sg[c] / (1-so[c]);
00510 so[c] = 0;
00511 }
00512 }
00513
00514 auto ixw = sw < 0;
00515 for (int c = 0; c < nc; ++c) {
00516 if (ixw[c]) {
00517 so[c] = so[c] / (1-sw[c]);
00518 sg[c] = sg[c] / (1-sw[c]);
00519 sw[c] = 0;
00520 }
00521 }
00522
00523 auto ixs = ss < 0;
00524 for (int c = 0; c < nc; ++c) {
00525 if (ixs[c]) {
00526 ss[c] = 0;
00527 }
00528 }
00529
00530
00531
00532
00533
00534
00535 so = V::Constant(nc,1.0) - sw - sg - ss;
00536
00537
00538 const double drmaxrel = drMaxRel();
00539 V rs;
00540 if (has_disgas_) {
00541 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
00542 const V drs = Base::isRs_ * dxvar;
00543 const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
00544 rs = rs_old - drs_limited;
00545 rs = rs.max(zero);
00546 }
00547 V rv;
00548 if (has_vapoil_) {
00549 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
00550 const V drv = Base::isRv_ * dxvar;
00551 const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
00552 rv = rv_old - drv_limited;
00553 rv = rv.max(zero);
00554 }
00555
00556 sd_.soMax = fluid_.satOilMax();
00557
00558
00559 const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
00560 auto watOnly = sw > (1 - epsilon);
00561
00562
00563 std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
00564 std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
00565
00566 if (has_disgas_) {
00567 const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
00568 const V rsSat = fluidRsSat(p, so, cells_);
00569 sd_.rsSat = ADB::constant(rsSat);
00570
00571 auto hasGas = (sg > 0 && Base::isRs_ == 0);
00572
00573
00574 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
00575 auto gasVaporized = ( (rs > rsSat * (1+epsilon) && Base::isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
00576 auto useSg = watOnly || hasGas || gasVaporized;
00577 for (int c = 0; c < nc; ++c) {
00578 if (useSg[c]) {
00579 rs[c] = rsSat[c];
00580 if (watOnly[c]) {
00581 so[c] = 0;
00582 sg[c] = 0;
00583 ss[c] = 0;
00584 rs[c] = 0;
00585 }
00586
00587 } else {
00588 hydroCarbonState[c] = HydroCarbonState::OilOnly;
00589 }
00590 }
00591 rs = rs.min(rsSat);
00592 }
00593
00594
00595 if (has_vapoil_) {
00596
00597
00598 const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
00599 const V gaspress = computeGasPressure(p, sw, so, sg);
00600 const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
00601 const V rvSat = fluidRvSat(gaspress, so, cells_);
00602 sd_.rvSat = ADB::constant(rvSat);
00603
00604
00605 auto hasOil = (so > 0 && Base::isRv_ == 0);
00606
00607
00608 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
00609 auto oilCondensed = ( (rv > rvSat * (1+epsilon) && Base::isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
00610 auto useSg = watOnly || hasOil || oilCondensed;
00611 for (int c = 0; c < nc; ++c) {
00612 if (useSg[c]) {
00613 rv[c] = rvSat[c];
00614 if (watOnly[c]) {
00615 so[c] = 0;
00616 sg[c] = 0;
00617 ss[c] = 0;
00618 rv[c] = 0;
00619 }
00620
00621 } else {
00622 hydroCarbonState[c] = HydroCarbonState::GasOnly;
00623 }
00624 }
00625 rv = rv.min(rvSat);
00626
00627 }
00628
00629
00630 if (has_solvent_) {
00631 auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
00632 std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
00633 }
00634
00635 for (int c = 0; c < nc; ++c) {
00636 reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
00637 }
00638
00639 for (int c = 0; c < nc; ++c) {
00640 reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
00641 }
00642
00643 if (active_[ Oil ]) {
00644 const int pos = pu.phase_pos[ Oil ];
00645 for (int c = 0; c < nc; ++c) {
00646 reservoir_state.saturation()[c*np + pos] = so[c];
00647 }
00648 }
00649
00650
00651 if (has_disgas_) {
00652 std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
00653 }
00654
00655 if (has_vapoil_) {
00656 std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
00657 }
00658
00659 wellModel().updateWellState(dwells, Base::dbhpMaxRel(), well_state);
00660
00661 for( auto w = 0; w < wells().number_of_wells; ++w ) {
00662 if (wells().type[w] == INJECTOR) {
00663 continue;
00664 }
00665
00666 for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf ) {
00667 int wc = wells().well_cells[perf];
00668 if ( (ss[wc] + sg[wc]) > 0) {
00669 well_state.solventFraction()[perf] = ss[wc] / (ss[wc] + sg[wc]);
00670 } else {
00671 well_state.solventFraction()[perf] = 0.0;
00672 }
00673 }
00674 }
00675
00676
00677
00678
00679 updatePhaseCondFromPrimalVariable(reservoir_state);
00680 }
00681
00682
00683
00684 template <class Grid>
00685 void
00686 BlackoilSolventModel<Grid>::computeMassFlux(const int actph ,
00687 const V& transi,
00688 const ADB& kr ,
00689 const ADB& mu ,
00690 const ADB& rho ,
00691 const ADB& phasePressure,
00692 const SolutionState& state)
00693 {
00694
00695 const int canonicalPhaseIdx = canph_[ actph ];
00696
00697 ADB kr_mod = kr;
00698 if (canonicalPhaseIdx == Gas) {
00699 if (has_solvent_) {
00700 const int nc = Opm::UgGridHelpers::numCells(grid_);
00701 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00702 const ADB zero = ADB::constant(V::Zero(nc));
00703 const V ones = V::Constant(nc, 1.0);
00704
00705 const ADB& ss = state.solvent_saturation;
00706 const ADB& sg = (active_[ Gas ]
00707 ? state.saturation[ pu.phase_pos[ Gas ] ]
00708 : zero);
00709 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00710 const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg));
00711
00712
00713 const std::vector<PhasePresence>& cond = phaseCondition();
00714 ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
00715 ADB rho_s = fluidDensity(Solvent,sd_.rq[solvent_pos_].b, state.rs, state.rv);
00716
00717
00718 ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;
00719 Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state);
00720
00721
00722 kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod;
00723
00724 }
00725 }
00726
00727 Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state);
00728
00729 }
00730
00731 template <class Grid>
00732 ADB
00733 BlackoilSolventModel<Grid>::fluidViscosity(const int phase,
00734 const ADB& p ,
00735 const ADB& temp ,
00736 const ADB& rs ,
00737 const ADB& rv ,
00738 const std::vector<PhasePresence>& cond) const
00739 {
00740 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00741 if (phase == Solvent) {
00742 if (!is_miscible_) {
00743 return solvent_props_.muSolvent(p, cells_);
00744 } else {
00745 return mu_eff_[solvent_pos_];
00746 }
00747
00748 } else {
00749 if (!is_miscible_) {
00750 return Base::fluidViscosity(phase, p, temp, rs, rv, cond);
00751 } else {
00752 return mu_eff_[pu.phase_pos[ phase ]];
00753 }
00754 }
00755 }
00756
00757 template <class Grid>
00758 ADB
00759 BlackoilSolventModel<Grid>::fluidReciprocFVF(const int phase,
00760 const ADB& p ,
00761 const ADB& temp ,
00762 const ADB& rs ,
00763 const ADB& rv ,
00764 const std::vector<PhasePresence>& cond) const
00765 {
00766 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00767 if (phase == Solvent) {
00768 if (!is_miscible_) {
00769 return solvent_props_.bSolvent(p, cells_);
00770 } else {
00771 return b_eff_[solvent_pos_];
00772 }
00773
00774 } else {
00775 if (!is_miscible_) {
00776 return Base::fluidReciprocFVF(phase, p, temp, rs, rv, cond);
00777 } else {
00778 return b_eff_[pu.phase_pos[ phase ]];
00779 }
00780 }
00781 }
00782
00783 template <class Grid>
00784 ADB
00785 BlackoilSolventModel<Grid>::fluidDensity(const int phase,
00786 const ADB& b,
00787 const ADB& rs,
00788 const ADB& rv) const
00789 {
00790 if (phase == Solvent && has_solvent_) {
00791 return solvent_props_.solventSurfaceDensity(cells_) * b;
00792 } else {
00793 return Base::fluidDensity(phase, b, rs, rv);
00794 }
00795 }
00796
00797 template <class Grid>
00798 std::vector<ADB>
00799 BlackoilSolventModel<Grid>::computeRelPerm(const SolutionState& state) const
00800 {
00801 using namespace Opm::AutoDiffGrid;
00802 const int nc = numCells(grid_);
00803
00804 const ADB zero = ADB::constant(V::Zero(nc));
00805
00806 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00807 const ADB& sw = (active_[ Water ]
00808 ? state.saturation[ pu.phase_pos[ Water ] ]
00809 : zero);
00810
00811 const ADB& so = (active_[ Oil ]
00812 ? state.saturation[ pu.phase_pos[ Oil ] ]
00813 : zero);
00814
00815 const ADB& sg = (active_[ Gas ]
00816 ? state.saturation[ pu.phase_pos[ Gas ] ]
00817 : zero);
00818
00819 if (has_solvent_) {
00820 const ADB& ss = state.solvent_saturation;
00821 if (is_miscible_) {
00822
00823 assert(active_[ Oil ]);
00824 assert(active_[ Gas ]);
00825
00826 std::vector<ADB> relperm = fluid_.relperm(sw, so, sg+ss, cells_);
00827
00828 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
00829 ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
00830 const ADB& po = state.canonical_phase_pressures[ Oil ];
00831 const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_)
00832 * solvent_props_.pressureMiscibilityFunction(po, cells_);
00833
00834 const ADB sn = ss + so + sg;
00835
00836
00837 const V sgcr = fluid_.scaledCriticalGasSaturations(cells_);
00838 const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_);
00839 const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
00840 const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
00841
00842 const V ones = V::Constant(nc, 1.0);
00843 ADB sor = misc * sorwmis + (ones - misc) * sogcr;
00844 ADB sgc = misc * sgcwmis + (ones - misc) * sgcr;
00845
00846 ADB ssg = ss + sg - sgc;
00847 const ADB sn_eff = sn - sor - sgc;
00848
00849
00850 Selector<double> negSsg_selector(ssg.value(), Selector<double>::LessZero);
00851 ssg = negSsg_selector.select(zero, ssg);
00852
00853
00854 Selector<double> zeroSn_selector(sn_eff.value(), Selector<double>::LessEqualZero);
00855 const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff);
00856
00857 const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
00858 const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
00859
00860 const V eps = V::Constant(nc, 1e-7);
00861 Selector<double> noOil_selector(so.value()-eps, Selector<double>::LessEqualZero);
00862 relperm[Gas] = (ones - misc) * relperm[Gas] + misc * mkrgt;
00863 relperm[Oil] = noOil_selector.select(relperm[Oil], (ones - misc) * relperm[Oil] + misc * mkro);
00864 return relperm;
00865 } else {
00866 return fluid_.relperm(sw, so, sg+ss, cells_);
00867 }
00868 } else {
00869 return fluid_.relperm(sw, so, sg, cells_);
00870 }
00871
00872 }
00873
00874
00875 template <class Grid>
00876 void
00877 BlackoilSolventModel<Grid>::computeEffectiveProperties(const SolutionState& state)
00878 {
00879
00880 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
00881 const int np = fluid_.numPhases();
00882 const int nc = Opm::UgGridHelpers::numCells(grid_);
00883 const ADB zero = ADB::constant(V::Zero(nc));
00884
00885 const ADB& pw = state.canonical_phase_pressures[pu.phase_pos[Water]];
00886 const ADB& po = state.canonical_phase_pressures[pu.phase_pos[Oil]];
00887 const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
00888 const std::vector<PhasePresence>& cond = phaseCondition();
00889
00890 const ADB mu_w = fluid_.muWat(pw, state.temperature, cells_);
00891 const ADB mu_o = fluid_.muOil(po, state.temperature, state.rs, cond, cells_);
00892 const ADB mu_g = fluid_.muGas(pg, state.temperature, state.rv, cond, cells_);
00893 const ADB mu_s = solvent_props_.muSolvent(pg,cells_);
00894 std::vector<ADB> viscosity(np + 1, ADB::null());
00895 viscosity[pu.phase_pos[Oil]] = mu_o;
00896 viscosity[pu.phase_pos[Gas]] = mu_g;
00897 viscosity[pu.phase_pos[Water]] = mu_w;
00898 viscosity[solvent_pos_] = mu_s;
00899
00900
00901 const ADB bw = fluid_.bWat(pw, state.temperature, cells_);
00902 const ADB bo = fluid_.bOil(po, state.temperature, state.rs, cond, cells_);
00903 const ADB bg = fluid_.bGas(pg, state.temperature, state.rv, cond, cells_);
00904 const ADB bs = solvent_props_.bSolvent(pg, cells_);
00905
00906 std::vector<ADB> density(np + 1, ADB::null());
00907 density[pu.phase_pos[Oil]] = fluidDensity(Oil, bo, state.rs, state.rv);
00908 density[pu.phase_pos[Gas]] = fluidDensity(Gas, bg, state.rs, state.rv);
00909 density[pu.phase_pos[Water]] = fluidDensity(Water, bw, state.rs, state.rv);
00910 density[solvent_pos_] = fluidDensity(Solvent, bs, state.rs, state.rv);
00911
00912
00913 const ADB& ss = state.solvent_saturation;
00914 const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ];
00915 const ADB& sg = (active_[ Gas ]
00916 ? state.saturation[ pu.phase_pos[ Gas ] ]
00917 : zero);
00918 const ADB& sw = (active_[ Water ]
00919 ? state.saturation[ pu.phase_pos[ Water ] ]
00920 : zero);
00921
00922 const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
00923 const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
00924
00925 std::vector<ADB> effective_saturations (np + 1, ADB::null());
00926 effective_saturations[pu.phase_pos[Oil]] = so - sorwmis;
00927 effective_saturations[pu.phase_pos[Gas]] = sg - sgcwmis;
00928 effective_saturations[pu.phase_pos[Water]] = sw;
00929 effective_saturations[solvent_pos_] = ss - sgcwmis;
00930
00931
00932 computeToddLongstaffMixing(viscosity, density, effective_saturations, po, pu);
00933
00934
00935 const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs);
00936 const ADB b_eff_g = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv);
00937 const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_);
00938
00939
00940 const V ones = V::Constant(nc, 1.0);
00941 const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
00942
00943 b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo;
00944 b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg;
00945 b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs;
00946
00947
00948 mu_eff_[pu.phase_pos[ Oil ]] = b_eff_[pu.phase_pos[ Oil ]] / (pmisc * b_eff_o / viscosity[pu.phase_pos[ Oil ]] + (ones - pmisc) * bo / mu_o);
00949 mu_eff_[pu.phase_pos[ Gas ]] = b_eff_[pu.phase_pos[ Gas ]] / (pmisc * b_eff_g / viscosity[pu.phase_pos[ Gas ]] + (ones - pmisc) * bg / mu_g);
00950 mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s);
00951
00952
00953 mu_eff_[pu.phase_pos[ Water ]] = mu_w;
00954 b_eff_[pu.phase_pos[ Water ]] = bw;
00955 }
00956
00957 template <class Grid>
00958 void
00959 BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const ADB po, const Opm::PhaseUsage pu)
00960 {
00961 const int nc = cells_.size();
00962 const V ones = V::Constant(nc, 1.0);
00963 const ADB zero = ADB::constant(V::Zero(nc));
00964
00965
00966 ADB so_eff = saturations[pu.phase_pos[ Oil ]];
00967 ADB sg_eff = saturations[pu.phase_pos[ Gas ]];
00968 ADB ss_eff = saturations[solvent_pos_];
00969
00970
00971 Selector<double> negative_selectorSo(so_eff.value(), Selector<double>::LessZero);
00972 Selector<double> negative_selectorSg(sg_eff.value(), Selector<double>::LessZero);
00973 Selector<double> negative_selectorSs(ss_eff.value(), Selector<double>::LessZero);
00974 so_eff = negative_selectorSo.select(zero, so_eff);
00975 sg_eff = negative_selectorSg.select(zero, sg_eff);
00976 ss_eff = negative_selectorSs.select(zero, ss_eff);
00977
00978
00979 const ADB mu_o = viscosity[pu.phase_pos[ Oil ]];
00980 const ADB mu_g = viscosity[pu.phase_pos[ Gas ]];
00981 const ADB mu_s = viscosity[solvent_pos_];
00982
00983 const ADB sn_eff = so_eff + sg_eff + ss_eff;
00984 const ADB sos_eff = so_eff + ss_eff;
00985 const ADB ssg_eff = ss_eff + sg_eff;
00986
00987
00988 Selector<double> zero_selectorSos(sos_eff.value(), Selector<double>::Zero);
00989 Selector<double> zero_selectorSsg(ssg_eff.value(), Selector<double>::Zero);
00990 Selector<double> zero_selectorSn(sn_eff.value(), Selector<double>::Zero);
00991
00992 const ADB mu_s_pow = pow(mu_s, 0.25);
00993 const ADB mu_o_pow = pow(mu_o, 0.25);
00994 const ADB mu_g_pow = pow(mu_g, 0.25);
00995
00996 const ADB mu_mos = zero_selectorSos.select(mu_o, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0));
00997 const ADB mu_msg = zero_selectorSsg.select(mu_g, mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0));
00998 const ADB mu_m = zero_selectorSn.select(mu_s, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow * mu_g_pow)
00999 + ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
01000
01001
01002
01003 const ADB mix_param_mu = solvent_props_.mixingParameterViscosity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
01004
01005 Selector<double> zero_selectorSs(ss_eff.value(), Selector<double>::Zero);
01006
01007 viscosity[pu.phase_pos[ Oil ]] = zero_selectorSs.select(mu_o, pow(mu_o,ones - mix_param_mu) * pow(mu_mos, mix_param_mu));
01008 viscosity[pu.phase_pos[ Gas ]] = zero_selectorSs.select(mu_g, pow(mu_g,ones - mix_param_mu) * pow(mu_msg, mix_param_mu));
01009 viscosity[solvent_pos_] = zero_selectorSs.select(mu_s, pow(mu_s,ones - mix_param_mu) * pow(mu_m, mix_param_mu));
01010
01011
01012 ADB& rho_o = density[pu.phase_pos[ Oil ]];
01013 ADB& rho_g = density[pu.phase_pos[ Gas ]];
01014 ADB& rho_s = density[solvent_pos_];
01015
01016
01017 const ADB mix_param_rho = solvent_props_.mixingParameterDensity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
01018
01019
01020
01021 const ADB mu_o_eff = pow(mu_o,ones - mix_param_rho) * pow(mu_mos, mix_param_rho);
01022 const ADB mu_g_eff = pow(mu_g,ones - mix_param_rho) * pow(mu_msg, mix_param_rho);
01023 const ADB mu_s_eff = pow(mu_s,ones - mix_param_rho) * pow(mu_m, mix_param_rho);
01024
01025 const ADB sog_eff = so_eff + sg_eff;
01026
01027 Selector<double> zero_selectorSog_eff(sog_eff.value(), Selector<double>::Zero);
01028 const ADB sof = zero_selectorSog_eff.select(zero , so_eff / sog_eff);
01029 const ADB sgf = zero_selectorSog_eff.select(zero , sg_eff / sog_eff);
01030
01031
01032 const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) );
01033 const ADB mu_o_eff_pow = pow(mu_o_eff, 0.25);
01034 const ADB mu_g_eff_pow = pow(mu_g_eff, 0.25);
01035 const ADB mu_s_eff_pow = pow(mu_s_eff, 0.25);
01036 const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow));
01037 const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow));
01038 const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow));
01039 const ADB rho_o_eff = (rho_o * sfraction_oe) + (rho_s * (ones - sfraction_oe));
01040 const ADB rho_g_eff = (rho_g * sfraction_ge) + (rho_s * (ones - sfraction_ge));
01041 const ADB rho_s_eff = (rho_s * sfraction_se) + (rho_g * sgf * (ones - sfraction_se)) + (rho_o * sof * (ones - sfraction_se));
01042
01043
01044
01045 Selector<double> unitGasSolventMobilityRatio_selector(mu_s.value() - mu_g.value(), Selector<double>::Zero);
01046 Selector<double> unitOilSolventMobilityRatio_selector(mu_s.value() - mu_o.value(), Selector<double>::Zero);
01047
01048
01049 const ADB rho_m = zero_selectorSn.select(zero, (rho_o * so_eff / sn_eff) + (rho_g * sg_eff / sn_eff) + (rho_s * ss_eff / sn_eff));
01050 const ADB rho_o_eff_simple = ((ones - mix_param_rho) * rho_o) + (mix_param_rho * rho_m);
01051 const ADB rho_g_eff_simple = ((ones - mix_param_rho) * rho_g) + (mix_param_rho * rho_m);
01052
01053
01054
01055 rho_o = zero_selectorSs.select(rho_o, unitOilSolventMobilityRatio_selector.select(rho_o_eff_simple, rho_o_eff) );
01056 rho_g = zero_selectorSs.select(rho_g, unitGasSolventMobilityRatio_selector.select(rho_g_eff_simple, rho_g_eff) );
01057 rho_s = zero_selectorSs.select(rho_s, rho_s_eff);
01058
01059
01060 }
01061
01062
01063 template <class Grid>
01064 std::vector<ADB>
01065 BlackoilSolventModel<Grid>::
01066 computePressures(const ADB& po,
01067 const ADB& sw,
01068 const ADB& so,
01069 const ADB& sg,
01070 const ADB& ss) const
01071 {
01072 std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
01073
01074 if (has_solvent_) {
01075
01076 std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
01077
01078
01079 const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
01080
01081
01082 const int nc = cells_.size();
01083 const V ones = V::Constant(nc, 1.0);
01084 pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
01085 }
01086 return pressures;
01087 }
01088
01089
01090
01091
01092 template <class Grid>
01093 std::vector<std::vector<double> >
01094 BlackoilSolventModel<Grid>::
01095 computeFluidInPlace(const ReservoirState& x,
01096 const std::vector<int>& fipnum)
01097 {
01098 if (has_solvent_ && is_miscible_ && b_eff_[0].size() == 0) {
01099
01100 WellState xw, xwdummy;
01101 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
01102 xw.init(&wells(), x, xwdummy, pu);
01103 SolutionState solstate = variableState(x, xw);
01104 computeEffectiveProperties(solstate);
01105 }
01106 return Base::computeFluidInPlace(x, fipnum);
01107 }
01108
01109
01110
01111 }
01112
01113
01114 #endif // OPM_BLACKOILSOLVENT_IMPL_HEADER_INCLUDED