fix bugs: polymer head should be dp*trans.

This commit is contained in:
Liu Ming 2014-10-20 13:44:12 +08:00
parent 89120ed57f
commit 88a1ec6e9f

View File

@ -1689,6 +1689,21 @@ namespace {
const ADB tr_mult = transMult(state.pressure); const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_); const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu; rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
ADB& head = rq_[phase].head;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
ADB dp = ops_.ngrad * phasePressure[canonicalPhaseIdx] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(dp);
}
head = transi*dp;
if (canonicalPhaseIdx == Water) { if (canonicalPhaseIdx == Water) {
if(has_polymer_) { if(has_polymer_) {
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
@ -1706,21 +1721,6 @@ namespace {
rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head; rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head;
} }
} }
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
ADB& head = rq_[phase].head;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
ADB dp = ops_.ngrad * phasePressure[canonicalPhaseIdx] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(dp);
}
head = transi*dp;
//head = transi*(ops_.ngrad * phasePressure) + gflux; //head = transi*(ops_.ngrad * phasePressure) + gflux;
UpwindSelector<double> upwind(grid_, ops_, head.value()); UpwindSelector<double> upwind(grid_, ops_, head.value());