Use pore volume to do the average when hydrocarbon is missing.

This commit is contained in:
Kai Bao
2021-03-05 11:24:35 +01:00
committed by Atgeirr Flø Rasmussen
parent d693c3da42
commit f9ead85435

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