diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp index dc44eabbb..3d77709ce 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp @@ -296,10 +296,12 @@ namespace Opm { computeRelPerm(const SolutionState& state) const; void - computeMassFlux(const V& transi, - const std::vector& kr , - const std::vector& phasePressure, + computeMassFlux(const int actph , + const V& transi, + const ADB& kr , + const ADB& p , const SolutionState& state ); + void computeCmax(PolymerBlackoilState& state); diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp index 624d13148..d2a986db4 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp @@ -872,9 +872,9 @@ namespace detail { // for each active phase. const V transi = subset(geo_.transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - computeMassFlux(transi, kr, state.canonical_phase_pressures, state); + // computeMassFlux(transi, kr, state.canonical_phase_pressures, state); for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { - // computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); + computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); residual_.material_balance_eq[ phaseIdx ] = pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + ops_.div*rq_[phaseIdx].mflux; @@ -1813,68 +1813,69 @@ namespace detail { return cp[Gas].value() + po; } - template void - FullyImplicitBlackoilPolymerSolver::computeMassFlux(const V& transi, - const std::vector& kr , - const std::vector& phasePressure, + FullyImplicitBlackoilPolymerSolver::computeMassFlux(const int actph , + const V& transi, + const ADB& kr , + const ADB& phasePressure, const SolutionState& state) { - for (int phase = 0; phase < fluid_.numPhases(); ++phase) { - const int canonicalPhaseIdx = canph_[phase]; - const std::vector cond = phaseCondition(); + const int canonicalPhaseIdx = canph_[ actph ]; - const ADB tr_mult = transMult(state.pressure); - const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_); + const std::vector cond = phaseCondition(); - rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu; + const ADB tr_mult = transMult(state.pressure); + const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_); - const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_); + rq_[ actph ].mob = tr_mult * kr / mu; - ADB& head = rq_[phase].head; + const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_); - // compute gravity potensial using the face average as in eclipse and MRST - const ADB rhoavg = ops_.caver * rho; + ADB& head = rq_[ actph ].head; - ADB dp = ops_.ngrad * phasePressure[canonicalPhaseIdx] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix())); + // compute gravity potensial using the face average as in eclipse and MRST + const ADB rhoavg = ops_.caver * rho; - if (use_threshold_pressure_) { - applyThresholdPressures(dp); + ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix())); + + if (use_threshold_pressure_) { + applyThresholdPressures(dp); + } + + head = transi*dp; + + if (canonicalPhaseIdx == Water) { + if(has_polymer_) { + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB mc = computeMc(state); + ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, + cmax, + kr, + state.saturation[canonicalPhaseIdx]); + ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data()); + rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc; + rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc; + rq_[poly_pos_].b = rq_[actph].b; + rq_[poly_pos_].head = rq_[actph].head; + UpwindSelector upwind(grid_, ops_, rq_[poly_pos_].head.value()); + rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head; } + } - head = transi*dp; - if (canonicalPhaseIdx == Water) { - if(has_polymer_) { - const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); - const ADB mc = computeMc(state); - ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, - cmax, - kr[canonicalPhaseIdx], - state.saturation[canonicalPhaseIdx]); - ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data()); - rq_[phase].mob = tr_mult * krw_eff * inv_wat_eff_visc; - rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc; - rq_[poly_pos_].b = rq_[phase].b; - rq_[poly_pos_].head = rq_[phase].head; - UpwindSelector upwind(grid_, ops_, rq_[poly_pos_].head.value()); - rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head; - } - } - //head = transi*(ops_.ngrad * phasePressure) + gflux; + //head = transi*(ops_.ngrad * phasePressure) + gflux; - UpwindSelector upwind(grid_, ops_, head.value()); + UpwindSelector upwind(grid_, ops_, head.value()); - const ADB& b = rq_[phase].b; - const ADB& mob = rq_[phase].mob; - rq_[phase].mflux = upwind.select(b * mob) * head; - } + const ADB& b = rq_[ actph ].b; + const ADB& mob = rq_[ actph ].mob; + rq_[ actph ].mflux = upwind.select(b * mob) * head; + // DUMP(rq_[ actph ].mob); + // DUMP(rq_[ actph ].mflux); } - - template void FullyImplicitBlackoilPolymerSolver::applyThresholdPressures(ADB& dp)