From f9ead854353b127f065cc743aba45458eec09a7a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 5 Mar 2021 11:24:35 +0100 Subject: [PATCH] Use pore volume to do the average when hydrocarbon is missing. --- opm/simulators/wells/RateConverter.hpp | 100 ++++++++++++++++--------- 1 file changed, 65 insertions(+), 35 deletions(-) diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 744e84b36..c22654691 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -446,6 +446,17 @@ namespace Opm { } + // quantities for pore volume average + std::unordered_map attributes_pv; + + // quantities for hydrocarbon volume average + std::unordered_map 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; + } } }