OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

refactoring computeMassFlux().

Browse Source 
To match the opm-autodiff implementation that only handles one phase
each time.

Polymer equation is handled with the water phase together, since it
effects the water phase.
This commit is contained in:
Kai Bao 2015-05-08 10:07:49 +02:00
parent 3bd0ff1885
commit b46dcc3b76
2 changed files with 51 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
View File

@ -296,10 +296,12 @@ namespace Opm {
computeRelPerm(const SolutionState& state) const; computeRelPerm(const SolutionState& state) const;
void void
computeMassFlux(const V& transi, computeMassFlux(const int actph ,
const std::vector<ADB>& kr , const V& transi,
const std::vector<ADB>& phasePressure, const ADB& kr ,
const ADB& p ,
const SolutionState& state ); const SolutionState& state );
void void
computeCmax(PolymerBlackoilState& state); computeCmax(PolymerBlackoilState& state);

95
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
View File

@ -872,9 +872,9 @@ namespace detail {
// for each active phase. // for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces); const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state); const std::vector<ADB> 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) { 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 ] = residual_.material_balance_eq[ phaseIdx ] =
pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux; + ops_.div*rq_[phaseIdx].mflux;
@ -1813,68 +1813,69 @@ namespace detail {
return cp[Gas].value() + po; return cp[Gas].value() + po;
} }
template<class T> template<class T>
void void
FullyImplicitBlackoilPolymerSolver<T>::computeMassFlux(const V& transi, FullyImplicitBlackoilPolymerSolver<T>::computeMassFlux(const int actph ,
const std::vector<ADB>& kr , const V& transi,
const std::vector<ADB>& phasePressure, const ADB& kr ,
const ADB& phasePressure,
const SolutionState& state) const SolutionState& state)
{ {
for (int phase = 0; phase < fluid_.numPhases(); ++phase) { const int canonicalPhaseIdx = canph_[ actph ];
const int canonicalPhaseIdx = canph_[phase];
const std::vector<PhasePresence> cond = phaseCondition();
const ADB tr_mult = transMult(state.pressure); const std::vector<PhasePresence> cond = phaseCondition();
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_);
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 ADB& head = rq_[ actph ].head;
const ADB rhoavg = ops_.caver * rho;
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_) { ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
applyThresholdPressures(dp);
}
head = transi*dp; if (use_threshold_pressure_) {
if (canonicalPhaseIdx == Water) { applyThresholdPressures(dp);
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<double> 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;
UpwindSelector<double> 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;
} }
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<double> 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;
UpwindSelector<double> upwind(grid_, ops_, head.value());
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<class T> template<class T>
void void
FullyImplicitBlackoilPolymerSolver<T>::applyThresholdPressures(ADB& dp) FullyImplicitBlackoilPolymerSolver<T>::applyThresholdPressures(ADB& dp)