Merge pull request #3099 from atgeirr/revive-using_porevolume_average

Fix hydrocarbon volume issue with no oil or gas
This commit is contained in:
Kai Bao 2021-03-09 12:15:06 +01:00 committed by GitHub
commit 0ad848658e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23

View File

@ -446,6 +446,17 @@ namespace Opm {
}
// quantities for pore volume average
std::unordered_map<RegionId, Attributes> attributes_pv;
// quantities for hydrocarbon volume average
std::unordered_map<RegionId, Attributes> attributes_hpv;
for (const auto& reg : rmap_.activeRegions()) {
attributes_pv.insert({reg, Attributes()});
attributes_hpv.insert({reg, Attributes()});
}
ElementContext elemCtx( simulator );
const auto& gridView = simulator.gridView();
const auto& comm = gridView.comm();
@ -477,49 +488,68 @@ namespace Opm {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
int reg = rmap_.region(cellIdx);
const int reg = rmap_.region(cellIdx);
assert(reg >= 0);
auto& ra = attr_.attributes(reg);
auto& p = ra.pressure;
auto& T = ra.temperature;
auto& rs = ra.rs;
auto& rv = ra.rv;
auto& pv = ra.pv;
auto& saltConcentration = ra.saltConcentration;
// sum p, rs, rv, and T.
double hydrocarbonPV = pv_cell*hydrocarbon;
if (hydrocarbonPV > 0) {
pv += hydrocarbonPV;
p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
rs += fs.Rs().value()*hydrocarbonPV;
rv += fs.Rv().value()*hydrocarbonPV;
T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
saltConcentration += fs.saltConcentration().value()*hydrocarbonPV;
const double hydrocarbonPV = pv_cell*hydrocarbon;
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
attr.rs += fs.Rs().value() * hydrocarbonPV;
attr.rv += fs.Rv().value() * hydrocarbonPV;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV;
}
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
attr.rs += fs.Rs().value() * pv_cell;
attr.rv += fs.Rv().value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
attr.saltConcentration += fs.saltConcentration().value() * pv_cell;
}
}
for (const auto& reg : rmap_.activeRegions()) {
auto& ra = attr_.attributes(reg);
auto& p = ra.pressure;
auto& T = ra.temperature;
auto& rs = ra.rs;
auto& rv = ra.rv;
auto& pv = ra.pv;
auto& saltConcentration = ra.saltConcentration;
// communicate sums
p = comm.sum(p);
T = comm.sum(T);
rs = comm.sum(rs);
rv = comm.sum(rv);
pv = comm.sum(pv);
saltConcentration = comm.sum(saltConcentration);
// compute average
p /= pv;
T /= pv;
rs /= pv;
rv /= pv;
saltConcentration /=pv;
const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
// TODO: should we have some epsilon here instead of zero?
if (hpv_sum > 0.) {
const auto& attri_hpv = attributes_hpv[reg];
const double p_hpv_sum = comm.sum(attri_hpv.pressure);
const double T_hpv_sum = comm.sum(attri_hpv.temperature);
const double rs_hpv_sum = comm.sum(attri_hpv.rs);
const double rv_hpv_sum = comm.sum(attri_hpv.rv);
const double sc_hpv_sum = comm.sum(attri_hpv.saltConcentration);
ra.pressure = p_hpv_sum / hpv_sum;
ra.temperature = T_hpv_sum / hpv_sum;
ra.rs = rs_hpv_sum / hpv_sum;
ra.rv = rv_hpv_sum / hpv_sum;
ra.pv = hpv_sum;
ra.saltConcentration = sc_hpv_sum / hpv_sum;
} else {
// using the pore volume to do the averaging
const auto& attri_pv = attributes_pv[reg];
const double pv_sum = comm.sum(attri_pv.pv);
assert(pv_sum > 0.);
const double p_pv_sum = comm.sum(attri_pv.pressure);
const double T_pv_sum = comm.sum(attri_pv.temperature);
const double rs_pv_sum = comm.sum(attri_pv.rs);
const double rv_pv_sum = comm.sum(attri_pv.rv);
const double sc_pv_sum = comm.sum(attri_pv.saltConcentration);
ra.pressure = p_pv_sum / pv_sum;
ra.temperature = T_pv_sum / pv_sum;
ra.rs = rs_pv_sum / pv_sum;
ra.rv = rv_pv_sum / pv_sum;
ra.pv = pv_sum;
ra.saltConcentration = sc_pv_sum / pv_sum;
}
}
}