From 6d3da3d2e0590224b8508cd2026b09313e287d81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 21 Jun 2021 16:15:10 +0200 Subject: [PATCH] Report Pressure Dependent Pore Volume in PRT File This commit distinguishes the reference condition pore volume from the dynamic, pressure (and/or temperature) dependent pore volume value. Previously we would report the latter as the 'PORV' value in the "Field Totals" and "FIPNUM region" reports, but this commit switches to reporting the former instead-mostly for compatibility. We still report the dynamic pore volume value, but now we report this on a line of its own, before the table, using one of the forms Field total pressure dependent pore volume = 12345 RM3 FIPNUM report region 1 pressure dependent pore volume = 123 RM3 --- ebos/eclgenericoutputblackoilmodule.cc | 292 ++++++++++++--------- ebos/eclgenericoutputblackoilmodule.hh | 11 +- ebos/ecloutputblackoilmodule.hh | 11 +- opm/simulators/wells/BlackoilWellModel.hpp | 11 +- opm/simulators/wells/WellState.cpp | 2 +- 5 files changed, 192 insertions(+), 135 deletions(-) diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index ff386865e..05daf6220 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -38,6 +38,7 @@ #include #include +#include #include #include #include @@ -827,33 +828,41 @@ doAllocBuffers(unsigned bufferSize, } } - outputFipRestart_ = false; - computeFip_ = false; + this->outputFipRestart_ = false; + this->computeFip_ = false; // Fluid in place for (const auto& phase : Inplace::phases()) { if (!substep || summaryConfig_.require3DField(EclString(phase))) { if (rstKeywords["FIP"] > 0) { rstKeywords["FIP"] = 0; - outputFipRestart_ = true; + this->outputFipRestart_ = true; } - fip_[phase].resize(bufferSize, 0.0); - computeFip_ = true; + + this->fip_[phase].resize(bufferSize, 0.0); + this->computeFip_ = true; + } + else { + this->fip_[phase].clear(); } - else - fip_[phase].clear(); } - if (!substep || summaryConfig_.hasKeyword("FPR") || summaryConfig_.hasKeyword("FPRP") || !this->RPRNodes_.empty()) { - fip_[Inplace::Phase::PoreVolume].resize(bufferSize, 0.0); - hydrocarbonPoreVolume_.resize(bufferSize, 0.0); - pressureTimesPoreVolume_.resize(bufferSize, 0.0); - pressureTimesHydrocarbonVolume_.resize(bufferSize, 0.0); + if (!substep || + this->summaryConfig_.hasKeyword("FPR") || + this->summaryConfig_.hasKeyword("FPRP") || + !this->RPRNodes_.empty()) + { + this->fip_[Inplace::Phase::PoreVolume].resize(bufferSize, 0.0); + this->dynamicPoreVolume_.resize(bufferSize, 0.0); + this->hydrocarbonPoreVolume_.resize(bufferSize, 0.0); + this->pressureTimesPoreVolume_.resize(bufferSize, 0.0); + this->pressureTimesHydrocarbonVolume_.resize(bufferSize, 0.0); } else { - hydrocarbonPoreVolume_.clear(); - pressureTimesPoreVolume_.clear(); - pressureTimesHydrocarbonVolume_.clear(); + this->dynamicPoreVolume_.clear(); + this->hydrocarbonPoreVolume_.clear(); + this->pressureTimesPoreVolume_.clear(); + this->pressureTimesHydrocarbonVolume_.clear(); } // Well RFT data @@ -1074,32 +1083,25 @@ void EclGenericOutputBlackoilModule:: fipUnitConvert_(std::unordered_map& fip) const { const UnitSystem& units = eclState_.getUnits(); - if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { - fip[Inplace::Phase::WATER] = unit::convert::to(fip[Inplace::Phase::WATER], unit::stb); - fip[Inplace::Phase::OIL] = unit::convert::to(fip[Inplace::Phase::OIL], unit::stb); - fip[Inplace::Phase::OilInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::OilInLiquidPhase], unit::stb); - fip[Inplace::Phase::OilInGasPhase] = unit::convert::to(fip[Inplace::Phase::OilInGasPhase], unit::stb); - fip[Inplace::Phase::GAS] = unit::convert::to(fip[Inplace::Phase::GAS], 1000*unit::cubic(unit::feet)); - fip[Inplace::Phase::GasInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::GasInLiquidPhase], 1000*unit::cubic(unit::feet)); - fip[Inplace::Phase::GasInGasPhase] = unit::convert::to(fip[Inplace::Phase::GasInGasPhase], 1000*unit::cubic(unit::feet)); - fip[Inplace::Phase::PoreVolume] = unit::convert::to(fip[Inplace::Phase::PoreVolume], unit::stb); - } - else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_LAB) { - Scalar scc = unit::cubic(prefix::centi * unit::meter); //standard cubic cm. - fip[Inplace::Phase::WATER] = unit::convert::to(fip[Inplace::Phase::WATER], scc); - fip[Inplace::Phase::OIL] = unit::convert::to(fip[Inplace::Phase::OIL], scc); - fip[Inplace::Phase::OilInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::OilInLiquidPhase], scc); - fip[Inplace::Phase::OilInGasPhase] = unit::convert::to(fip[Inplace::Phase::OilInGasPhase], scc); - fip[Inplace::Phase::GAS] = unit::convert::to(fip[Inplace::Phase::GAS], scc); - fip[Inplace::Phase::GasInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::GasInLiquidPhase], scc); - fip[Inplace::Phase::GasInGasPhase] = unit::convert::to(fip[Inplace::Phase::GasInGasPhase], scc); - fip[Inplace::Phase::PoreVolume] = unit::convert::to(fip[Inplace::Phase::PoreVolume], scc); - } - else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { - // nothing to do - } - else { - throw std::runtime_error("Unsupported unit type for fluid in place output."); + using M = UnitSystem::measure; + + const auto unit_map = std::unordered_map { + {Inplace::Phase::WATER, M::liquid_surface_volume}, + {Inplace::Phase::OIL, M::liquid_surface_volume}, + {Inplace::Phase::OilInLiquidPhase, M::liquid_surface_volume}, + {Inplace::Phase::OilInGasPhase, M::liquid_surface_volume}, + {Inplace::Phase::GAS, M::gas_surface_volume}, + {Inplace::Phase::GasInLiquidPhase, M::gas_surface_volume}, + {Inplace::Phase::GasInGasPhase, M::gas_surface_volume}, + {Inplace::Phase::PoreVolume, M::volume}, + {Inplace::Phase::DynamicPoreVolume, M::volume}, + }; + + for (auto& [phase, value] : fip) { + auto unitPos = unit_map.find(phase); + if (unitPos != unit_map.end()) { + value = units.from_si(unitPos->second, value); + } } } @@ -1107,20 +1109,8 @@ template void EclGenericOutputBlackoilModule:: pressureUnitConvert_(Scalar& pav) const { - const UnitSystem& units = eclState_.getUnits(); - if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { - pav = unit::convert::to(pav, unit::psia); - } - else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { - pav = unit::convert::to(pav, unit::barsa); - } - else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_LAB) { - pav = unit::convert::to(pav, unit::atm); - - } - else { - throw std::runtime_error("Unsupported unit type for fluid in place output."); - } + pav = this->eclState_.getUnits() + .from_si(UnitSystem::measure::pressure, pav); } template @@ -1133,11 +1123,25 @@ outputRegionFluidInPlace_(std::unordered_map oip, return; // don't output FIPNUM report if the region has no porv. - if (cip[Inplace::Phase::PoreVolume] == 0) + if (! (cip[Inplace::Phase::PoreVolume] > Scalar{0})) return; const UnitSystem& units = eclState_.getUnits(); std::ostringstream ss; + + ss << '\n'; + if (reg == 0) { + ss << "Field total"; + } + else { + ss << "FIPNUM report region " << reg; + } + + ss << " pressure dependent pore volume = " + << std::fixed << std::setprecision(0) + << cip[Inplace::Phase::DynamicPoreVolume] << ' ' + << units.name(UnitSystem::measure::volume) << "\n\n"; + if (reg == 0) { ss << " ===================================================\n" << " : Field Totals :\n"; @@ -1304,12 +1308,12 @@ isOutputCreationDirective_(const std::string& keyword) template Scalar EclGenericOutputBlackoilModule:: pressureAverage_(const Scalar& pressurePvHydrocarbon, - const Scalar& pvHydrocarbon, - const Scalar& pressurePv, - const Scalar& pv, - bool hydrocarbon) + const Scalar& pvHydrocarbon, + const Scalar& pressurePv, + const Scalar& pv, + const bool hydrocarbon) { - if (pvHydrocarbon > 1e-10 && hydrocarbon) + if (hydrocarbon && (pvHydrocarbon > 1e-10)) return pressurePvHydrocarbon / pvHydrocarbon; return pressurePv / pv; @@ -1319,20 +1323,25 @@ template typename EclGenericOutputBlackoilModule::ScalarBuffer EclGenericOutputBlackoilModule:: pressureAverage_(const ScalarBuffer& pressurePvHydrocarbon, - const ScalarBuffer& pvHydrocarbon, - const ScalarBuffer& pressurePv, - const ScalarBuffer& pv, - bool hydrocarbon) + const ScalarBuffer& pvHydrocarbon, + const ScalarBuffer& pressurePv, + const ScalarBuffer& pv, + const bool hydrocarbon) { - size_t size = pressurePvHydrocarbon.size(); + const std::size_t size = pressurePvHydrocarbon.size(); assert(pvHydrocarbon.size() == size); assert(pressurePv.size() == size); assert(pv.size() == size); ScalarBuffer fraction(size, 0.0); - for (size_t i = 0; i < size; ++i) { - fraction[i] = pressureAverage_(pressurePvHydrocarbon[i], pvHydrocarbon[i], pressurePv[i], pv[i], hydrocarbon); + for (std::size_t i = 0; i < size; ++i) { + fraction[i] = pressureAverage_(pressurePvHydrocarbon[i], + pvHydrocarbon[i], + pressurePv[i], + pv[i], + hydrocarbon); } + return fraction; } @@ -1421,13 +1430,15 @@ outputFipLogImpl(const Inplace& inplace) const current_values[phase] = inplace.get(phase); } + current_values[Inplace::Phase::DynamicPoreVolume] = + inplace.get(Inplace::Phase::DynamicPoreVolume); fipUnitConvert_(initial_values); fipUnitConvert_(current_values); pressureUnitConvert_(fieldHydroCarbonPoreVolumeAveragedPressure); - outputRegionFluidInPlace_(initial_values, - current_values, + outputRegionFluidInPlace_(std::move(initial_values), + std::move(current_values), fieldHydroCarbonPoreVolumeAveragedPressure); } @@ -1439,6 +1450,10 @@ outputFipLogImpl(const Inplace& inplace) const initial_values[phase] = this->initialInplace_->get("FIPNUM", phase, reg); current_values[phase] = inplace.get("FIPNUM", phase, reg); } + + current_values[Inplace::Phase::DynamicPoreVolume] = + inplace.get("FIPNUM", Inplace::Phase::DynamicPoreVolume, reg); + fipUnitConvert_(initial_values); fipUnitConvert_(current_values); @@ -1449,7 +1464,9 @@ outputFipLogImpl(const Inplace& inplace) const inplace.get("FIPNUM", Inplace::Phase::PoreVolume, reg), true); pressureUnitConvert_(regHydroCarbonPoreVolumeAveragedPressure); - outputRegionFluidInPlace_(initial_values, current_values, regHydroCarbonPoreVolumeAveragedPressure, reg); + outputRegionFluidInPlace_(std::move(initial_values), + std::move(current_values), + regHydroCarbonPoreVolumeAveragedPressure, reg); } } @@ -1466,33 +1483,55 @@ template void EclGenericOutputBlackoilModule:: update(Inplace& inplace, const std::string& region_name, - Inplace::Phase phase, - std::size_t ntFip, - const std::vector& values) + const Inplace::Phase phase, + const std::size_t ntFip, + const ScalarBuffer& values) { - double sum = 0; - for (std::size_t region_number = 0; region_number < ntFip; region_number++) { - inplace.add( region_name, phase, region_number + 1, values[region_number] ); - sum += values[region_number]; + double sum = 0.0; + for (std::size_t region_number = 0; region_number < ntFip; ++region_number) { + const auto rval = static_cast(values[region_number]); + inplace.add(region_name, phase, region_number + 1, rval); + sum += rval; } - inplace.add( phase, sum ); + inplace.add(phase, sum); } template void EclGenericOutputBlackoilModule:: makeRegionSum(Inplace& inplace, const std::string& region_name, - const Comm& comm) + const Comm& comm) const { const auto& region = this->regions_.at(region_name); - std::size_t ntFip = this->regionMax(region, comm); + const std::size_t ntFip = this->regionMax(region, comm); - update(inplace, region_name, Inplace::Phase::PressurePV, ntFip, this->regionSum(this->pressureTimesPoreVolume_, region, ntFip, comm)); - update(inplace, region_name, Inplace::Phase::HydroCarbonPV, ntFip, this->regionSum(this->hydrocarbonPoreVolume_, region, ntFip, comm)); - update(inplace, region_name, Inplace::Phase::PressureHydroCarbonPV, ntFip, this->regionSum(this->pressureTimesHydrocarbonVolume_, region, ntFip, comm)); + auto update_inplace = + [&inplace, ®ion, ®ion_name, &comm, ntFip, this] + (const Inplace::Phase phase, + const std::vector& value) + { + update(inplace, region_name, phase, ntFip, + this->regionSum(value, region, ntFip, comm)); + }; - for (const auto& phase : Inplace::phases()) - update(inplace, region_name, phase, ntFip, this->regionSum(this->fip_[phase], region, ntFip, comm)); + update_inplace(Inplace::Phase::PressurePV, + this->pressureTimesPoreVolume_); + + update_inplace(Inplace::Phase::HydroCarbonPV, + this->hydrocarbonPoreVolume_); + + update_inplace(Inplace::Phase::PressureHydroCarbonPV, + this->pressureTimesHydrocarbonVolume_); + + update_inplace(Inplace::Phase::DynamicPoreVolume, + this->dynamicPoreVolume_); + + for (const auto& phase : Inplace::phases()) { + auto fipPos = this->fip_.find(phase); + if (fipPos != this->fip_.end()) { + update_inplace(phase, fipPos->second); + } + } } template @@ -1501,9 +1540,8 @@ accumulateRegionSums(const Comm& comm) { Inplace inplace; - for (const auto& [region_name, _] : this->regions_) { - (void)_; - makeRegionSum(inplace, region_name, comm); + for (const auto& region : this->regions_) { + makeRegionSum(inplace, region.first, comm); } // The first time the outputFipLog function is run we store the inplace values in @@ -1537,54 +1575,68 @@ updateSummaryRegionValues(const Inplace& inplace, // The field summary vectors should only use the FIPNUM based region sum. { for (const auto& phase : Inplace::phases()) { - std::string key = "F" + EclString(phase); - if (summaryConfig_.hasKeyword(key)) + const std::string key = "F" + EclString(phase); + if (this->summaryConfig_.hasKeyword(key)) { miscSummaryData[key] = inplace.get(phase); + } } - if (summaryConfig_.hasKeyword("FOE") && this->initialInplace_) + if (this->summaryConfig_.hasKeyword("FOE") && this->initialInplace_) { miscSummaryData["FOE"] = inplace.get(Inplace::Phase::OIL) / this->initialInplace_.value().get(Inplace::Phase::OIL); + } - if (summaryConfig_.hasKeyword("FPR")) - miscSummaryData["FPR"] = pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV), - inplace.get(Inplace::Phase::HydroCarbonPV), - inplace.get(Inplace::Phase::PressurePV), - inplace.get(Inplace::Phase::PoreVolume), - true); + if (this->summaryConfig_.hasKeyword("FPR")) { + miscSummaryData["FPR"] = + pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV), + inplace.get(Inplace::Phase::HydroCarbonPV), + inplace.get(Inplace::Phase::PressurePV), + inplace.get(Inplace::Phase::PoreVolume), + true); + } - - if (summaryConfig_.hasKeyword("FPRP")) - miscSummaryData["FPRP"] = pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV), - inplace.get(Inplace::Phase::HydroCarbonPV), - inplace.get(Inplace::Phase::PressurePV), - inplace.get(Inplace::Phase::PoreVolume), - false); + if (this->summaryConfig_.hasKeyword("FPRP")) { + miscSummaryData["FPRP"] = + pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV), + inplace.get(Inplace::Phase::HydroCarbonPV), + inplace.get(Inplace::Phase::PressurePV), + inplace.get(Inplace::Phase::PoreVolume), + false); + } } // The region summary vectors should loop through the FIPxxx regions to // support the RPR__xxx summary keywords. { + auto get_vector = [&inplace] + (const auto& node, + const Inplace::Phase phase) + { + return inplace.get_vector(node.fip_region(), phase); + }; + for (const auto& phase : Inplace::phases()) { for (const auto& node : this->regionNodes_.at(phase)) - regionData[node.keyword()] = inplace.get_vector(node.fip_region(), phase); + regionData[node.keyword()] = get_vector(node, phase); } - // The exact same quantity is calculated for RPR and RPRP - is that correct? - for (const auto& node : this->RPRNodes_) - regionData[node.keyword()] = pressureAverage_(inplace.get_vector(node.fip_region(), Inplace::Phase::PressureHydroCarbonPV), - inplace.get_vector(node.fip_region(), Inplace::Phase::HydroCarbonPV), - inplace.get_vector(node.fip_region(), Inplace::Phase::PressurePV), - inplace.get_vector(node.fip_region(), Inplace::Phase::PoreVolume), - true); + for (const auto& node : this->RPRNodes_) { + regionData[node.keyword()] = + pressureAverage_(get_vector(node, Inplace::Phase::PressureHydroCarbonPV), + get_vector(node, Inplace::Phase::HydroCarbonPV), + get_vector(node, Inplace::Phase::PressurePV), + get_vector(node, Inplace::Phase::PoreVolume), + true); + } - - for (const auto& node : this->RPRPNodes_) - regionData[node.keyword()] = pressureAverage_(inplace.get_vector(node.fip_region(), Inplace::Phase::PressureHydroCarbonPV), - inplace.get_vector(node.fip_region(), Inplace::Phase::HydroCarbonPV), - inplace.get_vector(node.fip_region(), Inplace::Phase::PressurePV), - inplace.get_vector(node.fip_region(), Inplace::Phase::PoreVolume), - false); + for (const auto& node : this->RPRPNodes_) { + regionData[node.keyword()] = + pressureAverage_(get_vector(node, Inplace::Phase::PressureHydroCarbonPV), + get_vector(node, Inplace::Phase::HydroCarbonPV), + get_vector(node, Inplace::Phase::PressurePV), + get_vector(node, Inplace::Phase::PoreVolume), + false); + } } } diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh index 64bd10ee9..d184f6bce 100644 --- a/ebos/eclgenericoutputblackoilmodule.hh +++ b/ebos/eclgenericoutputblackoilmodule.hh @@ -261,7 +261,7 @@ protected: void makeRegionSum(Inplace& inplace, const std::string& region_name, - const Comm& comm); + const Comm& comm) const; Inplace accumulateRegionSums(const Comm& comm); @@ -285,7 +285,7 @@ protected: // Sum Fip values over regions. static ScalarBuffer regionSum(const ScalarBuffer& property, const std::vector& regionId, - size_t maxNumberOfRegions, + const std::size_t maxNumberOfRegions, const Comm& comm); static int regionMax(const std::vector& region, @@ -293,9 +293,9 @@ protected: static void update(Inplace& inplace, const std::string& region_name, - Inplace::Phase phase, - std::size_t ntFip, - const std::vector& values); + const Inplace::Phase phase, + const std::size_t ntFip, + const ScalarBuffer& values); static Scalar sum(const ScalarBuffer& v); @@ -332,6 +332,7 @@ protected: ScalarBuffer hydrocarbonPoreVolume_; ScalarBuffer pressureTimesPoreVolume_; ScalarBuffer pressureTimesHydrocarbonVolume_; + ScalarBuffer dynamicPoreVolume_; ScalarBuffer oilPressure_; ScalarBuffer temperature_; ScalarBuffer rs_; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 3e3542697..6f3408f65 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -623,15 +623,18 @@ private: // PORV, MINPV and friends). Also note that because of this, the porosity // returned by the intensive quantities can be outside of the physical // range [0, 1] in pathetic cases. - const double pv = - elemCtx.simulator().model().dofTotalVolume(globalDofIdx) - * intQuants.porosity().value(); + const auto totVolume = + elemCtx.simulator().model().dofTotalVolume(globalDofIdx); + const double pv = totVolume * intQuants.porosity().value(); if (!this->pressureTimesHydrocarbonVolume_.empty() && !this->pressureTimesPoreVolume_.empty()) { assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size()); assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size()); - this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = pv; + this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = + totVolume * intQuants.referencePorosity(); + + this->dynamicPoreVolume_[globalDofIdx] = pv; Scalar hydrocarbon = 0.0; if (FluidSystem::phaseIsActive(oilPhaseIdx)) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 04208ddea..c92e05787 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -222,11 +222,12 @@ namespace Opm { data::Wells wellData() const { - auto wsrpt = this->wellState().report(UgGridHelpers::globalCell(grid()), - [this](const int well_ndex) -> bool - { - return this->wasDynamicallyShutThisTimeStep(well_ndex); - }); + auto wsrpt = this->wellState() + .report(UgGridHelpers::globalCell(grid()), + [this](const int well_ndex) -> bool + { + return this->wasDynamicallyShutThisTimeStep(well_ndex); + }); this->assignWellGuideRates(wsrpt); this->assignShutConnections(wsrpt, this->reportStepIndex()); diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 6a44ad9be..798195fcf 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -492,7 +492,7 @@ void WellState::gatherVectorsOnRoot(const std::vector& from_co data::Wells WellState::report(const int* globalCellIdxMap, - const std::function& wasDynamicallyClosed) const + const std::function& wasDynamicallyClosed) const { if (this->numWells() == 0) return {};