Added missing code for polymer model for output

This commit is contained in:
babrodtk 2016-09-15 11:40:36 +02:00
parent bbd2575e00
commit 537c5d71b8
2 changed files with 14 additions and 9 deletions

View File

@ -257,7 +257,7 @@ namespace Opm {
/// Computing the water velocity without shear-thinning for the cell faces.
/// The water velocity will be used for shear-thinning calculation.
void computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
void computeWaterShearVelocityFaces(const V& transi, const std::vector<typename Base::ReservoirResidualQuant>& rq,
const std::vector<ADB>& phasePressure, const SolutionState& state,
std::vector<double>& water_vel, std::vector<double>& visc_mult);

View File

@ -291,14 +291,19 @@ namespace Opm {
// Set up the common parts of the mass balance equations
// for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state);
{
const std::vector<ADB> kr = computeRelPerm(state);
for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
}
}
if (has_plyshlog_) {
std::vector<double> water_vel;
std::vector<double> visc_mult;
computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult);
computeWaterShearVelocityFaces(transi, sd_.rq, state.canonical_phase_pressures, state, water_vel, visc_mult);
if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
// std::cerr << " failed in calculating the shear-multiplier " << std::endl;
OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
@ -307,9 +312,9 @@ namespace Opm {
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
const std::vector<PhasePresence>& cond = phaseCondition();
const ADB mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
const ADB rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
@ -585,7 +590,7 @@ namespace Opm {
template<class Grid>
void
BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<typename Base::ReservoirResidualQuant>& rq,
const std::vector<ADB>& phasePressure, const SolutionState& state,
std::vector<double>& water_vel, std::vector<double>& visc_mult)
{
@ -598,7 +603,7 @@ namespace Opm {
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
sd_.rq[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
sd_.rq[phase].mob = tr_mult * rq[phase].kr / mu;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
@ -617,7 +622,7 @@ namespace Opm {
const ADB mc = computeMc(state);
ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
cmax,
kr[canonicalPhaseIdx]);
rq[phase].kr);
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;