20 #ifndef OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED 21 #define OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED 23 #include <opm/autodiff/BlackoilSolventModel.hpp> 25 #include <opm/autodiff/AutoDiffBlock.hpp> 26 #include <opm/autodiff/AutoDiffHelpers.hpp> 27 #include <opm/autodiff/GridHelpers.hpp> 28 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp> 29 #include <opm/autodiff/GeoProps.hpp> 30 #include <opm/autodiff/WellDensitySegmented.hpp> 32 #include <opm/core/grid.h> 33 #include <opm/core/linalg/LinearSolverInterface.hpp> 34 #include <opm/core/linalg/ParallelIstlInformation.hpp> 35 #include <opm/core/props/rock/RockCompressibility.hpp> 36 #include <opm/common/ErrorMacros.hpp> 37 #include <opm/common/Exceptions.hpp> 38 #include <opm/parser/eclipse/Units/Units.hpp> 39 #include <opm/core/well_controls.h> 40 #include <opm/core/utility/parameters/ParameterGroup.hpp> 55 int solventPos(
const PU& pu)
57 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
59 for (
int phase = 0; phase < maxnp; ++phase) {
60 if (pu.phase_used[phase]) {
77 const RockCompressibility* rock_comp_props,
81 std::shared_ptr< const EclipseState > eclState,
82 const bool has_disgas,
83 const bool has_vapoil,
84 const bool terminal_output,
85 const bool has_solvent,
86 const bool is_miscible)
87 :
Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
88 eclState, has_disgas, has_vapoil, terminal_output),
89 has_solvent_(has_solvent),
90 solvent_pos_(detail::solventPos(fluid.phaseUsage())),
91 solvent_props_(solvent_props),
92 is_miscible_(is_miscible)
100 Base::material_name_.push_back(
"Solvent");
101 assert(solvent_pos_ == fluid_.
numPhases());
102 residual_.matbalscale.resize(fluid_.
numPhases() + 1, 0.0031);
104 wellModel().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
116 template <
class Gr
id>
120 Base::makeConstantState(state);
121 state.solvent_saturation =
ADB::constant(state.solvent_saturation.value());
128 template <
class Gr
id>
130 BlackoilSolventModel<Grid>::variableStateInitials(
const ReservoirState& x,
131 const WellState& xw)
const 133 std::vector<V> vars0 = Base::variableStateInitials(x, xw);
134 assert(
int(vars0.size()) == fluid_.numPhases() + 2);
138 const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL );
139 const int nc = solvent_saturation.size();
140 const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
144 auto solvent_pos = vars0.begin() + fluid_.numPhases();
145 assert (not solvent_saturation.empty());
146 assert(solvent_pos == vars0.end() - 2);
147 vars0.insert(solvent_pos, ss);
156 template <
class Gr
id>
158 BlackoilSolventModel<Grid>::variableStateIndices()
const 160 std::vector<int> ind = Base::variableStateIndices();
161 assert(ind.size() == 5);
165 ind[Solvent] = fluid_.numPhases();
176 template <
class Gr
id>
177 typename BlackoilSolventModel<Grid>::SolutionState
178 BlackoilSolventModel<Grid>::variableStateExtractVars(
const ReservoirState& x,
179 const std::vector<int>& indices,
180 std::vector<ADB>& vars)
const 186 const int nc = Opm::AutoDiffGrid::numCells(grid_);
187 const Opm::PhaseUsage pu = fluid_.phaseUsage();
189 SolutionState state(fluid_.numPhases());
192 state.pressure = std::move(vars[indices[Pressure]]);
195 const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
202 if (active_[ Water ]) {
203 state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
204 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
208 state.solvent_saturation = std::move(vars[indices[Solvent]]);
209 so -= state.solvent_saturation;
212 if (active_[ Gas ]) {
215 const ADB& xvar = vars[indices[Xvar]];
216 ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
217 sg = Base::isSg_*xvar + Base::isRv_*so;
220 if (active_[ Oil ]) {
222 state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation);
223 sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
225 state.rs = (1-Base::isRs_)*sd_.rsSat + Base::isRs_*xvar;
227 state.rs = sd_.rsSat;
229 sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
231 state.rv = (1-Base::isRv_)*sd_.rvSat + Base::isRv_*xvar;
233 state.rv = sd_.rvSat;
238 if (active_[ Oil ]) {
240 state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
244 wellModel().variableStateExtractWellsVars(indices, vars, state);
253 template <
class Gr
id>
255 BlackoilSolventModel<Grid>::computeAccum(
const SolutionState& state,
260 computeEffectiveProperties(state);
262 Base::computeAccum(state, aix);
266 const ADB& press = state.pressure;
267 const ADB& ss = state.solvent_saturation;
268 const ADB pv_mult = poroMult(press);
269 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
271 const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
272 const std::vector<PhasePresence>& cond = phaseCondition();
273 sd_.rq[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
274 sd_.rq[solvent_pos_].accum[aix] = pv_mult * sd_.rq[solvent_pos_].b * ss;
282 template <
class Gr
id>
284 BlackoilSolventModel<Grid>::
285 assembleMassBalanceEq(
const SolutionState& state)
288 Base::assembleMassBalanceEq(state);
291 residual_.material_balance_eq[ solvent_pos_ ] =
292 pvdt_ * (sd_.rq[solvent_pos_].accum[1] - sd_.rq[solvent_pos_].accum[0])
293 + ops_.div*sd_.rq[solvent_pos_].mflux;
298 template <
class Gr
id>
300 BlackoilSolventModel<Grid>::updateEquationsScaling()
302 Base::updateEquationsScaling();
303 assert(MaxNumPhases + 1 == residual_.matbalscale.size());
305 const ADB& temp_b = sd_.rq[solvent_pos_].b;
306 ADB::V B = 1. / temp_b.value();
308 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
310 const ParallelISTLInformation& real_info =
311 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
312 double B_global_sum = 0;
313 real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
314 residual_.matbalscale[solvent_pos_] = B_global_sum / Base::global_nc_;
319 residual_.matbalscale[solvent_pos_] = B.mean();
324 template <
class Gr
id>
325 void BlackoilSolventModel<Grid>::addWellContributionToMassBalanceEq(
const std::vector<ADB>& cq_s,
326 const SolutionState& state,
333 Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
336 const int nperf = wells().well_connpos[wells().number_of_wells];
337 const int nc = Opm::AutoDiffGrid::numCells(grid_);
339 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
341 const ADB& ss = state.solvent_saturation;
342 const ADB& sg = (active_[ Gas ]
343 ? state.saturation[ pu.phase_pos[ Gas ] ]
346 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
347 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
348 ADB F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
350 const int nw = wells().number_of_wells;
351 V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
353 V isProducer = V::Zero(nperf);
354 V ones = V::Constant(nperf,1.0);
355 for (
int w = 0; w < nw; ++w) {
356 if(wells().type[w] == PRODUCER) {
357 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
358 isProducer[perf] = 1;
363 const ADB& rs_perfcells =
subset(state.rs, well_cells);
364 const ADB& rv_perfcells =
subset(state.rv, well_cells);
365 int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
366 int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
369 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);
373 residual_.material_balance_eq[solvent_pos_] -=
superset(cq_s_solvent, well_cells, nc);
377 residual_.material_balance_eq[gas_pos] +=
superset(cq_s_solvent, well_cells, nc);
386 template <
class Gr
id>
390 ReservoirState& reservoir_state,
391 WellState& well_state)
399 const int np = fluid_.numPhases();
400 const int nc = numCells(grid_);
402 assert(null.size() == 0);
403 const V zero = V::Zero(nc);
404 const V ones = V::Constant(nc,1.0);
409 const V dsw = active_[Water] ?
subset(dx,
Span(nc, 1, varstart)) : null;
410 varstart += dsw.size();
412 const V dxvar = active_[Gas] ?
subset(dx,
Span(nc, 1, varstart)): null;
413 varstart += dxvar.size();
415 const V dss = has_solvent_ ?
subset(dx,
Span(nc, 1, varstart)) : null;
416 varstart += dss.size();
419 const V dwells =
subset(dx,
Span(wellModel().numWellVars(), 1, varstart));
420 varstart += dwells.size();
422 assert(varstart == dx.size());
425 const double dpmaxrel = dpMaxRel();
426 const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
427 const V absdpmax = dpmaxrel*p_old.abs();
428 const V dp_limited =
sign(dp) * dp.abs().min(absdpmax);
429 const V p = (p_old - dp_limited).max(zero);
430 std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
434 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
435 const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
436 const double dsmax = dsMax();
451 maxVal = dsw.abs().max(maxVal);
457 dsg = Base::isSg_ * dxvar - Base::isRv_ * dsw;
458 maxVal = dsg.abs().max(maxVal);
462 maxVal = dss.abs().max(maxVal);
468 maxVal = dso.abs().max(maxVal);
470 V step = dsmax/maxVal;
473 if (active_[Water]) {
474 const int pos = pu.phase_pos[ Water ];
475 const V sw_old = s_old.col(pos);
476 sw = sw_old - step * dsw;
480 const int pos = pu.phase_pos[ Gas ];
481 const V sg_old = s_old.col(pos);
482 sg = sg_old - step * dsg;
485 auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
486 const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
487 ss = ss_old - step * dss;
490 const int pos = pu.phase_pos[ Oil ];
491 const V so_old = s_old.col(pos);
492 so = so_old - step * dso;
497 for (
int c = 0; c < nc; ++c) {
499 sw[c] = sw[c] / (1-sg[c]);
500 so[c] = so[c] / (1-sg[c]);
506 for (
int c = 0; c < nc; ++c) {
508 sw[c] = sw[c] / (1-so[c]);
509 sg[c] = sg[c] / (1-so[c]);
515 for (
int c = 0; c < nc; ++c) {
517 so[c] = so[c] / (1-sw[c]);
518 sg[c] = sg[c] / (1-sw[c]);
524 for (
int c = 0; c < nc; ++c) {
535 so = V::Constant(nc,1.0) - sw - sg - ss;
538 const double drmaxrel = drMaxRel();
541 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
542 const V drs = Base::isRs_ * dxvar;
543 const V drs_limited =
sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
544 rs = rs_old - drs_limited;
549 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
550 const V drv = Base::isRv_ * dxvar;
551 const V drv_limited =
sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
552 rv = rv_old - drv_limited;
556 sd_.soMax = fluid_.satOilMax();
559 const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
560 auto watOnly = sw > (1 - epsilon);
563 std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
564 std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
567 const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
568 const V rsSat = fluidRsSat(p, so, cells_);
571 auto hasGas = (sg > 0 && Base::isRs_ == 0);
574 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
575 auto gasVaporized = ( (rs > rsSat * (1+epsilon) && Base::isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
576 auto useSg = watOnly || hasGas || gasVaporized;
577 for (
int c = 0; c < nc; ++c) {
588 hydroCarbonState[c] = HydroCarbonState::OilOnly;
598 const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
599 const V gaspress = computeGasPressure(p, sw, so, sg);
600 const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
601 const V rvSat = fluidRvSat(gaspress, so, cells_);
605 auto hasOil = (so > 0 && Base::isRv_ == 0);
608 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
609 auto oilCondensed = ( (rv > rvSat * (1+epsilon) && Base::isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
610 auto useSg = watOnly || hasOil || oilCondensed;
611 for (
int c = 0; c < nc; ++c) {
622 hydroCarbonState[c] = HydroCarbonState::GasOnly;
631 auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
632 std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
635 for (
int c = 0; c < nc; ++c) {
636 reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
639 for (
int c = 0; c < nc; ++c) {
640 reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
643 if (active_[ Oil ]) {
644 const int pos = pu.phase_pos[ Oil ];
645 for (
int c = 0; c < nc; ++c) {
646 reservoir_state.saturation()[c*np + pos] = so[c];
652 std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
656 std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
659 wellModel().updateWellState(dwells, Base::dbhpMaxRel(), well_state);
661 for(
auto w = 0; w < wells().number_of_wells; ++w ) {
662 if (wells().type[w] == INJECTOR) {
666 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf ) {
667 int wc = wells().well_cells[perf];
668 if ( (ss[wc] + sg[wc]) > 0) {
669 well_state.solventFraction()[perf] = ss[wc] / (ss[wc] + sg[wc]);
671 well_state.solventFraction()[perf] = 0.0;
679 updatePhaseCondFromPrimalVariable(reservoir_state);
684 template <
class Gr
id>
691 const ADB& phasePressure,
692 const SolutionState& state)
695 const int canonicalPhaseIdx = canph_[ actph ];
698 if (canonicalPhaseIdx == Gas) {
700 const int nc = Opm::UgGridHelpers::numCells(grid_);
701 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
703 const V ones = V::Constant(nc, 1.0);
705 const ADB& ss = state.solvent_saturation;
706 const ADB& sg = (active_[ Gas ]
707 ? state.saturation[ pu.phase_pos[ Gas ] ]
710 const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg));
713 const std::vector<PhasePresence>& cond = phaseCondition();
714 ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
715 ADB rho_s = fluidDensity(Solvent,sd_.rq[solvent_pos_].b, state.rs, state.rv);
718 ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;
719 Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state);
722 kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod;
727 Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state);
731 template <
class Gr
id>
733 BlackoilSolventModel<Grid>::fluidViscosity(
const int phase,
738 const std::vector<PhasePresence>& cond)
const 740 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
741 if (phase == Solvent) {
743 return solvent_props_.muSolvent(p, cells_);
745 return mu_eff_[solvent_pos_];
750 return Base::fluidViscosity(phase, p, temp, rs, rv, cond);
752 return mu_eff_[pu.phase_pos[ phase ]];
757 template <
class Gr
id>
759 BlackoilSolventModel<Grid>::fluidReciprocFVF(
const int phase,
764 const std::vector<PhasePresence>& cond)
const 766 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
767 if (phase == Solvent) {
769 return solvent_props_.bSolvent(p, cells_);
771 return b_eff_[solvent_pos_];
776 return Base::fluidReciprocFVF(phase, p, temp, rs, rv, cond);
778 return b_eff_[pu.phase_pos[ phase ]];
783 template <
class Gr
id>
785 BlackoilSolventModel<Grid>::fluidDensity(
const int phase,
790 if (phase == Solvent && has_solvent_) {
791 return solvent_props_.solventSurfaceDensity(cells_) * b;
793 return Base::fluidDensity(phase, b, rs, rv);
797 template <
class Gr
id>
799 BlackoilSolventModel<Grid>::computeRelPerm(
const SolutionState& state)
const 802 const int nc = numCells(grid_);
806 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
807 const ADB& sw = (active_[ Water ]
808 ? state.saturation[ pu.phase_pos[ Water ] ]
811 const ADB& so = (active_[ Oil ]
812 ? state.saturation[ pu.phase_pos[ Oil ] ]
815 const ADB& sg = (active_[ Gas ]
816 ? state.saturation[ pu.phase_pos[ Gas ] ]
820 const ADB& ss = state.solvent_saturation;
823 assert(active_[ Oil ]);
824 assert(active_[ Gas ]);
826 std::vector<ADB> relperm = fluid_.relperm(sw, so, sg+ss, cells_);
828 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
829 ADB F_solvent = zero_selector.
select(ss, ss / (ss + sg));
830 const ADB& po = state.canonical_phase_pressures[ Oil ];
831 const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_)
832 * solvent_props_.pressureMiscibilityFunction(po, cells_);
834 const ADB sn = ss + so + sg;
837 const V sgcr = fluid_.scaledCriticalGasSaturations(cells_);
838 const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_);
839 const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
840 const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
842 const V ones = V::Constant(nc, 1.0);
843 ADB sor = misc * sorwmis + (ones - misc) * sogcr;
844 ADB sgc = misc * sgcwmis + (ones - misc) * sgcr;
846 ADB ssg = ss + sg - sgc;
847 const ADB sn_eff = sn - sor - sgc;
850 Selector<double> negSsg_selector(ssg.value(), Selector<double>::LessZero);
851 ssg = negSsg_selector.select(zero, ssg);
854 Selector<double> zeroSn_selector(sn_eff.value(), Selector<double>::LessEqualZero);
855 const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff);
857 const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
858 const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
860 const V eps = V::Constant(nc, 1e-7);
861 Selector<double> noOil_selector(so.value()-eps, Selector<double>::LessEqualZero);
862 relperm[Gas] = (ones - misc) * relperm[Gas] + misc * mkrgt;
863 relperm[Oil] = noOil_selector.select(relperm[Oil], (ones - misc) * relperm[Oil] + misc * mkro);
866 return fluid_.relperm(sw, so, sg+ss, cells_);
869 return fluid_.relperm(sw, so, sg, cells_);
875 template <
class Gr
id>
877 BlackoilSolventModel<Grid>::computeEffectiveProperties(
const SolutionState& state)
880 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
881 const int np = fluid_.numPhases();
882 const int nc = Opm::UgGridHelpers::numCells(grid_);
885 const ADB& pw = state.canonical_phase_pressures[pu.phase_pos[Water]];
886 const ADB& po = state.canonical_phase_pressures[pu.phase_pos[Oil]];
887 const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
888 const std::vector<PhasePresence>& cond = phaseCondition();
890 const ADB mu_w = fluid_.muWat(pw, state.temperature, cells_);
891 const ADB mu_o = fluid_.muOil(po, state.temperature, state.rs, cond, cells_);
892 const ADB mu_g = fluid_.muGas(pg, state.temperature, state.rv, cond, cells_);
893 const ADB mu_s = solvent_props_.muSolvent(pg,cells_);
894 std::vector<ADB> viscosity(np + 1,
ADB::null());
895 viscosity[pu.phase_pos[Oil]] = mu_o;
896 viscosity[pu.phase_pos[Gas]] = mu_g;
897 viscosity[pu.phase_pos[Water]] = mu_w;
898 viscosity[solvent_pos_] = mu_s;
901 const ADB bw = fluid_.bWat(pw, state.temperature, cells_);
902 const ADB bo = fluid_.bOil(po, state.temperature, state.rs, cond, cells_);
903 const ADB bg = fluid_.bGas(pg, state.temperature, state.rv, cond, cells_);
904 const ADB bs = solvent_props_.bSolvent(pg, cells_);
906 std::vector<ADB> density(np + 1,
ADB::null());
907 density[pu.phase_pos[Oil]] = fluidDensity(Oil, bo, state.rs, state.rv);
908 density[pu.phase_pos[Gas]] = fluidDensity(Gas, bg, state.rs, state.rv);
909 density[pu.phase_pos[Water]] = fluidDensity(Water, bw, state.rs, state.rv);
910 density[solvent_pos_] = fluidDensity(Solvent, bs, state.rs, state.rv);
913 const ADB& ss = state.solvent_saturation;
914 const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ];
915 const ADB& sg = (active_[ Gas ]
916 ? state.saturation[ pu.phase_pos[ Gas ] ]
918 const ADB& sw = (active_[ Water ]
919 ? state.saturation[ pu.phase_pos[ Water ] ]
922 const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
923 const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
925 std::vector<ADB> effective_saturations (np + 1,
ADB::null());
926 effective_saturations[pu.phase_pos[Oil]] = so - sorwmis;
927 effective_saturations[pu.phase_pos[Gas]] = sg - sgcwmis;
928 effective_saturations[pu.phase_pos[Water]] = sw;
929 effective_saturations[solvent_pos_] = ss - sgcwmis;
932 computeToddLongstaffMixing(viscosity, density, effective_saturations, po, pu);
935 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);
936 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);
937 const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_);
940 const V ones = V::Constant(nc, 1.0);
941 const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
943 b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo;
944 b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg;
945 b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs;
948 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);
949 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);
950 mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s);
953 mu_eff_[pu.phase_pos[ Water ]] = mu_w;
954 b_eff_[pu.phase_pos[ Water ]] = bw;
957 template <
class Gr
id>
959 BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density,
const std::vector<ADB>& saturations,
const ADB po,
const Opm::PhaseUsage pu)
961 const int nc = cells_.
size();
962 const V ones = V::Constant(nc, 1.0);
966 ADB so_eff = saturations[pu.phase_pos[ Oil ]];
967 ADB sg_eff = saturations[pu.phase_pos[ Gas ]];
968 ADB ss_eff = saturations[solvent_pos_];
971 Selector<double> negative_selectorSo(so_eff.value(), Selector<double>::LessZero);
972 Selector<double> negative_selectorSg(sg_eff.value(), Selector<double>::LessZero);
973 Selector<double> negative_selectorSs(ss_eff.value(), Selector<double>::LessZero);
974 so_eff = negative_selectorSo.select(zero, so_eff);
975 sg_eff = negative_selectorSg.select(zero, sg_eff);
976 ss_eff = negative_selectorSs.select(zero, ss_eff);
979 const ADB mu_o = viscosity[pu.phase_pos[ Oil ]];
980 const ADB mu_g = viscosity[pu.phase_pos[ Gas ]];
981 const ADB mu_s = viscosity[solvent_pos_];
983 const ADB sn_eff = so_eff + sg_eff + ss_eff;
984 const ADB sos_eff = so_eff + ss_eff;
985 const ADB ssg_eff = ss_eff + sg_eff;
988 Selector<double> zero_selectorSos(sos_eff.value(), Selector<double>::Zero);
989 Selector<double> zero_selectorSsg(ssg_eff.value(), Selector<double>::Zero);
990 Selector<double> zero_selectorSn(sn_eff.value(), Selector<double>::Zero);
992 const ADB mu_s_pow =
pow(mu_s, 0.25);
993 const ADB mu_o_pow =
pow(mu_o, 0.25);
994 const ADB mu_g_pow =
pow(mu_g, 0.25);
996 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));
997 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));
998 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)
999 + ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
1003 const ADB mix_param_mu = solvent_props_.mixingParameterViscosity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
1005 Selector<double> zero_selectorSs(ss_eff.value(), Selector<double>::Zero);
1007 viscosity[pu.phase_pos[ Oil ]] = zero_selectorSs.select(mu_o,
pow(mu_o,ones - mix_param_mu) *
pow(mu_mos, mix_param_mu));
1008 viscosity[pu.phase_pos[ Gas ]] = zero_selectorSs.select(mu_g,
pow(mu_g,ones - mix_param_mu) *
pow(mu_msg, mix_param_mu));
1009 viscosity[solvent_pos_] = zero_selectorSs.select(mu_s,
pow(mu_s,ones - mix_param_mu) *
pow(mu_m, mix_param_mu));
1012 ADB& rho_o = density[pu.phase_pos[ Oil ]];
1013 ADB& rho_g = density[pu.phase_pos[ Gas ]];
1014 ADB& rho_s = density[solvent_pos_];
1017 const ADB mix_param_rho = solvent_props_.mixingParameterDensity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
1021 const ADB mu_o_eff =
pow(mu_o,ones - mix_param_rho) *
pow(mu_mos, mix_param_rho);
1022 const ADB mu_g_eff =
pow(mu_g,ones - mix_param_rho) *
pow(mu_msg, mix_param_rho);
1023 const ADB mu_s_eff =
pow(mu_s,ones - mix_param_rho) *
pow(mu_m, mix_param_rho);
1025 const ADB sog_eff = so_eff + sg_eff;
1027 Selector<double> zero_selectorSog_eff(sog_eff.value(), Selector<double>::Zero);
1028 const ADB sof = zero_selectorSog_eff.select(zero , so_eff / sog_eff);
1029 const ADB sgf = zero_selectorSog_eff.select(zero , sg_eff / sog_eff);
1032 const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) );
1033 const ADB mu_o_eff_pow =
pow(mu_o_eff, 0.25);
1034 const ADB mu_g_eff_pow =
pow(mu_g_eff, 0.25);
1035 const ADB mu_s_eff_pow =
pow(mu_s_eff, 0.25);
1036 const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow));
1037 const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow));
1038 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));
1039 const ADB rho_o_eff = (rho_o * sfraction_oe) + (rho_s * (ones - sfraction_oe));
1040 const ADB rho_g_eff = (rho_g * sfraction_ge) + (rho_s * (ones - sfraction_ge));
1041 const ADB rho_s_eff = (rho_s * sfraction_se) + (rho_g * sgf * (ones - sfraction_se)) + (rho_o * sof * (ones - sfraction_se));
1045 Selector<double> unitGasSolventMobilityRatio_selector(mu_s.value() - mu_g.value(), Selector<double>::Zero);
1046 Selector<double> unitOilSolventMobilityRatio_selector(mu_s.value() - mu_o.value(), Selector<double>::Zero);
1049 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));
1050 const ADB rho_o_eff_simple = ((ones - mix_param_rho) * rho_o) + (mix_param_rho * rho_m);
1051 const ADB rho_g_eff_simple = ((ones - mix_param_rho) * rho_g) + (mix_param_rho * rho_m);
1055 rho_o = zero_selectorSs.select(rho_o, unitOilSolventMobilityRatio_selector.select(rho_o_eff_simple, rho_o_eff) );
1056 rho_g = zero_selectorSs.select(rho_g, unitGasSolventMobilityRatio_selector.select(rho_g_eff_simple, rho_g_eff) );
1057 rho_s = zero_selectorSs.select(rho_s, rho_s_eff);
1063 template <
class Gr
id>
1065 BlackoilSolventModel<Grid>::
1066 computePressures(
const ADB& po,
1070 const ADB& ss)
const 1072 std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
1076 std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
1079 const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
1082 const int nc = cells_.
size();
1083 const V ones = V::Constant(nc, 1.0);
1084 pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
1092 template <
class Gr
id>
1093 std::vector<std::vector<double> >
1094 BlackoilSolventModel<Grid>::
1095 computeFluidInPlace(
const ReservoirState& x,
1096 const std::vector<int>& fipnum)
1098 if (has_solvent_ && is_miscible_ && b_eff_[0].size() == 0) {
1100 WellState xw, xwdummy;
1101 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1102 xw.init(&wells(), x, xwdummy, pu);
1103 SolutionState solstate = variableState(x, xw);
1104 computeEffectiveProperties(solstate);
1106 return Base::computeFluidInPlace(x, fipnum);
1114 #endif // OPM_BLACKOILSOLVENT_IMPL_HEADER_INCLUDED static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
BlackoilSolventModel(const typename Base::ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const SolventPropsAdFromDeck &solvent_props, const StandardWellsSolvent &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output, const bool has_solvent, const bool is_miscible)
Construct the model.
Definition: BlackoilSolventModel_impl.hpp:73
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
Class for handling the standard well model for solvent model.
Definition: StandardWellsSolvent.hpp:32
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilSolventModel_impl.hpp:389
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:364
AutoDiffBlock< Scalar > pow(const AutoDiffBlock< Scalar > &base, const double exponent)
Computes the value of base raised to the power of exponent.
Definition: AutoDiffBlock.hpp:636
std::vector< ADB > material_balance_eq
The material_balance_eq vector has one element for each active phase, each of which has size equal to...
Definition: LinearisedBlackoilResidual.hpp:54
Definition: GridHelpers.cpp:27
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
Definition: AutoDiffHelpers.hpp:623
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
Definition: SolventPropsAdFromDeck.hpp:37
ADB select(const ADB &x1, const ADB &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:415
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
A model implementation for three-phase black oil with one extra component.
Definition: BlackoilSolventModel.hpp:38
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719