Merge pull request #808 from qilicun/fix-pav

Correct the formulation for PAV calculations.
This commit is contained in:
Atgeirr Flø Rasmussen 2016-09-07 08:52:18 +02:00 committed by GitHub
commit f78155fe23
4 changed files with 71 additions and 36 deletions

View File

@ -2310,9 +2310,9 @@ namespace detail {
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np); const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
state.pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1)); state.pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
state.temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1)); state.temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
state.saturation[Water] = ADB::constant(s.col(Water)); state.saturation[Water] = active_[Water] ? ADB::constant(s.col(Water)) : ADB::null();
state.saturation[Oil] = ADB::constant(s.col(Oil)); state.saturation[Oil] = active_[Oil] ? ADB::constant(s.col(Oil)) : ADB::constant(V::Zero(nc));
state.saturation[Gas] = ADB::constant(s.col(Gas)); state.saturation[Gas] = active_[Gas] ? ADB::constant(s.col(Gas)) : ADB::constant(V::Zero(nc));
state.rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1)); state.rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
state.rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1)); state.rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1));
state.canonical_phase_pressures = computePressures(state.pressure, state.canonical_phase_pressures = computePressures(state.pressure,
@ -2324,8 +2324,6 @@ namespace detail {
const std::vector<PhasePresence> cond = phaseCondition(); const std::vector<PhasePresence> cond = phaseCondition();
const ADB pv_mult = poroMult(state.pressure); const ADB pv_mult = poroMult(state.pressure);
const V& pv = geo_.poreVolume(); const V& pv = geo_.poreVolume();
const int maxnp = Opm::BlackoilPhases::MaxNumPhases; const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
@ -2356,16 +2354,25 @@ namespace detail {
} }
} }
// compute PAV and PORV or every regions. // compute PAV and PORV for every regions.
const V hydrocarbon = state.saturation[Oil].value() + state.saturation[Gas].value();
V hcpv = V::Zero(nc);
V pres = V::Zero(nc);
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) { if (fipnum[c] != 0) {
hcpv[fipnum[c]-1] += pv[c] * hydrocarbon[c];
pres[fipnum[c]-1] += pv[c] * state.pressure.value()[c];
values[fipnum[c]-1][5] += pv[c]; values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * state.pressure.value()[c]; values[fipnum[c]-1][6] += pv[c] * state.pressure.value()[c] * hydrocarbon[c];
} }
} }
for (auto& val : values) { for (int reg = 0; reg < dims; ++reg) {
val[6] = val[6] / val[5]; if (hcpv[reg] != 0) {
values[reg][6] /= hcpv[reg];
} else {
values[reg][6] = pres[reg] / values[reg][5];
}
} }
return values; return values;

View File

@ -163,7 +163,7 @@ namespace Opm
std::vector<V>& fip); std::vector<V>& fip);
V V
FIPTotals(const std::vector<V>& fip, const std::vector<double>& press); FIPTotals(const std::vector<V>& fip, const ReservoirState& state);
void void
outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg); outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg);

View File

@ -271,8 +271,8 @@ namespace Opm
std::vector<V> COIP; std::vector<V> COIP;
COIP = solver->computeFluidInPlace(state, fipnum); COIP = solver->computeFluidInPlace(state, fipnum);
FIPUnitConvert(eclipse_state_->getUnits(), COIP); FIPUnitConvert(eclipse_state_->getUnits(), COIP);
V OOIP_totals = FIPTotals(OOIP, state.pressure()); V OOIP_totals = FIPTotals(OOIP, state);
V COIP_totals = FIPTotals(COIP, state.pressure()); V COIP_totals = FIPTotals(COIP, state);
outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0); outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) { for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1); outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
@ -682,7 +682,7 @@ namespace Opm
template <class Implementation> template <class Implementation>
V V
SimulatorBase<Implementation>::FIPTotals(const std::vector<V>& fip, const std::vector<double>& press) SimulatorBase<Implementation>::FIPTotals(const std::vector<V>& fip, const ReservoirState& state)
{ {
V totals(V::Zero(7)); V totals(V::Zero(7));
for (int i = 0; i < 5; ++i) { for (int i = 0; i < 5; ++i) {
@ -690,9 +690,16 @@ namespace Opm
totals[i] += fip[reg][i]; totals[i] += fip[reg][i];
} }
} }
const V p = Eigen::Map<const V>(& press[0], press.size()); const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = state.numPhases();
const PhaseUsage& pu = props_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
const V so = pu.phase_used[BlackoilPhases::Liquid] ? V(s.col(BlackoilPhases::Liquid)) : V::Zero(nc);
const V sg = pu.phase_used[BlackoilPhases::Vapour] ? V(s.col(BlackoilPhases::Vapour)) : V::Zero(nc);
const V hydrocarbon = so + sg;
const V p = Eigen::Map<const V>(& state.pressure()[0], nc);
totals[5] = geo_.poreVolume().sum(); totals[5] = geo_.poreVolume().sum();
totals[6] = unit::convert::to((p * geo_.poreVolume()).sum() / totals[5], unit::barsa); totals[6] = unit::convert::to((p * geo_.poreVolume() * hydrocarbon).sum() / ((geo_.poreVolume() * hydrocarbon).sum()), unit::barsa);
return totals; return totals;
} }
@ -705,28 +712,32 @@ namespace Opm
{ {
std::ostringstream ss; std::ostringstream ss;
if (!reg) { if (!reg) {
ss << " ==================================================\n" ss << " ===================================================\n"
<< " : Field Totals :\n"; << " : Field Totals :\n";
} else { } else {
ss << " ==================================================\n" ss << " ===================================================\n"
<< " : FIPNUM report region " << " : FIPNUM report region "
<< std::setw(2) << reg << " :\n"; << std::setw(2) << reg << " :\n";
} }
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n" ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
<< std::fixed << std::setprecision(0) << std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n" << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
<< " : Pressure is weighted by hydrocarbon pore voulme:\n" if (!reg) {
<< " : Porv volume are taken at reference conditions :\n" ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n"; << " : Porv volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
} }
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :" ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
<< std::fixed << std::setprecision(0) << std::fixed << std::setprecision(0)
<< " : PORV =" << std::setprecision(14) << cip[5] << " STB :" << " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
<< " : Pressure is weighted by hydrocarbon pore voulme :" if (!reg) {
<< " : Pore volume are taken at reference conditions :" ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
<< " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n"; << " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
} }
ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n" ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n" << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"

View File

@ -564,16 +564,33 @@ namespace {
const int dims = *std::max_element(fipnum.begin(), fipnum.end()); const int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<V> values(dims, V::Zero(5)); std::vector<V> values(dims, V::Zero(7));
V hcpv = V::Zero(nc);
for (int d = 0; d < dims; ++d) { V pres = V::Zero(nc);
for (int c = 0; c < nc; ++c) {
for (int i = 0; i < 5; ++i) { for (int i = 0; i < 5; ++i) {
if (fipnum[c] == d) { for (int c = 0; c < nc; ++c) {
values[d][i] += fip[c][i]; if (fipnum[c] != 0) {
values[fipnum[c]-1][i] += fip[i][c];
} }
} }
} }
// compute PAV and PORV or every regions.
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
hcpv[fipnum[c]-1] += pv[c] * s.col(Oil)[c];
pres[fipnum[c]-1] += pv[c] * state.pressure.value()[c];
values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * state.pressure.value()[c] * s.col(Oil)[c];
}
}
for (int reg = 0; reg < dims; ++reg) {
if (hcpv[reg] != 0) {
values[reg][6] /= hcpv[reg];
} else {
values[reg][6] = pres[reg] / values[reg][5];
}
} }
return values; return values;