diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 91dcb28ad..331389c03 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -27,6 +27,11 @@ #ifndef EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH #define EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH +#include +#include +#include +#include + #include #include @@ -65,6 +70,8 @@ struct ForceDisableFluidInPlaceOutput { namespace Opm { + + // forward declaration template class EcfvDiscretization; @@ -114,7 +121,22 @@ class EclOutputBlackOilModule GasInPlaceInGasPhase = 6, //GIPG PoreVolume = 7, //PV }; - static const int numFipValues = PoreVolume + 1 ; + static const int numFipTypes = PoreVolume + 1 ; + + + static std::string EclString(int fip_type) { + switch(static_cast(fip_type)) { + case FipDataType::WaterInPlace: return "WIP"; + case FipDataType::OilInPlace: return "OIP"; + case FipDataType::GasInPlace: return "GIP"; + case FipDataType::OilInPlaceInLiquidPhase: return "OIPL"; + case FipDataType::OilInPlaceInGasPhase: return "OIPG"; + case FipDataType::GasInPlaceInLiquidPhase: return "GIPL"; + case FipDataType::GasInPlaceInGasPhase: return "GIPG"; + case FipDataType::PoreVolume: return "PV"; + } + throw std::logic_error("fip_type: " + std::to_string(fip_type) + " not recognized"); + } }; struct WellProdDataType { @@ -181,17 +203,42 @@ class EclOutputBlackOilModule static const int numWCNames = 3; }; + struct RegionSum { + std::size_t ntFip; + std::array regFipValues; + ScalarBuffer regPressurePv; + ScalarBuffer regPvHydrocarbon; + ScalarBuffer regPressurePvHydrocarbon; + + std::array fieldFipValues; + Scalar fieldPressurePv; + Scalar fieldPvHydrocarbon; + Scalar fieldPressurePvHydrocarbon; + }; + public: template EclOutputBlackOilModule(const Simulator& simulator, const CollectDataToIORankType& collectToIORank) : simulator_(simulator) { - createLocalFipnum_(); - - // Summary output is for all steps const Opm::SummaryConfig summaryConfig = simulator_.vanguard().summaryConfig(); + const auto& fp = simulator_.vanguard().eclState().fieldProps(); + + this->regions_["FIPNUM"] = fp.get_int("FIPNUM"); + for (const auto& region : summaryConfig.fip_regions()) + this->regions_[region] = fp.get_int(region); + + for (auto& region_pair : this->regions_) + createLocalRegion_(region_pair.second); + + this->RPRNodes_ = summaryConfig.keywords("RPR*"); + this->RPRPNodes_ = summaryConfig.keywords("RPRP*"); + + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) { + std::string key_pattern = "R" + FipDataType::EclString(fip_type) + "*"; + this->regionNodes_[fip_type] = summaryConfig.keywords(key_pattern); + } - // Initialize block output for (const auto& node: summaryConfig) { if (node.category() == SummaryConfigNode::Category::Block) { if (collectToIORank.isCartIdxOnThisRank(node.number() - 1)) { @@ -244,8 +291,8 @@ public: computeFip_ = false; // Fluid in place - for (int i = 0; i 0) { rstKeywords["FIP"] = 0; outputFipRestart_ = true; @@ -256,7 +303,7 @@ public: else fip_[i].clear(); } - if (!substep || summaryConfig.hasKeyword("FPR") || summaryConfig.hasKeyword("FPRP") || summaryConfig.hasKeyword("RPR")) { + if (!substep || summaryConfig.hasKeyword("FPR") || summaryConfig.hasKeyword("FPRP") || !this->RPRNodes_.empty()) { fip_[FipDataType::PoreVolume].resize(bufferSize, 0.0); hydrocarbonPoreVolume_.resize(bufferSize, 0.0); pressureTimesPoreVolume_.resize(bufferSize, 0.0); @@ -1025,9 +1072,9 @@ public: // Fluid in place - for (int i = 0; i 0) { - sol.insert(fipEnumToString_(i), + sol.insert(FipDataType::EclString(i), Opm::UnitSystem::measure::volume, fip_[i], Opm::data::TargetType::SUMMARY); @@ -1044,106 +1091,204 @@ public: } } + int regionMax(const std::vector& region) { + const auto max_value = region.empty() ? 0 : *std::max_element(region.begin(), region.end()); + return this->simulator_.gridView().comm().max(max_value); + } + + + RegionSum makeRegionSum(const std::vector& region, bool is_fipnum) { + RegionSum rsum; + rsum.ntFip = this->regionMax(region); + + // sum values over each region + rsum.regPressurePv = this->regionSum(this->pressureTimesPoreVolume_, region, rsum.ntFip); + rsum.regPvHydrocarbon = this->regionSum(this->hydrocarbonPoreVolume_, region, rsum.ntFip); + rsum.regPressurePvHydrocarbon = this->regionSum(pressureTimesHydrocarbonVolume_, region, rsum.ntFip); + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) + rsum.regFipValues[fip_type] = this->regionSum(this->fip_[fip_type], region, rsum.ntFip); + + + + // sum all region values to compute the field total + rsum.fieldPressurePv = sum(rsum.regPressurePv); + rsum.fieldPvHydrocarbon = sum(rsum.regPvHydrocarbon); + rsum.fieldPressurePvHydrocarbon = sum(rsum.regPressurePvHydrocarbon); + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) { + const auto& regionFip = rsum.regFipValues[fip_type]; + rsum.fieldFipValues[fip_type] = sum(regionFip); + } + + if (is_fipnum) { + // The first time the outputFipLog function is run we store the inplace values in + // the initialInplace_ member. This has at least two problems: + // + // o We really want the *initial* value - now we get the value after + // the first timestep. + // + // o For restarted runs this is obviously wrong. + // + // Finally it is of course not desirable to mutate state in an output + // routine. + if (!this->regionInitialInplace_) { + this->regionInitialInplace_.emplace(); + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) + this->regionInitialInplace_.value()[fip_type] = rsum.regFipValues[fip_type]; + } + + if (!this->fieldInitialInplace_) + this->fieldInitialInplace_ = rsum.fieldFipValues; + } + return rsum; + } + + + std::unordered_map accumulateRegionSums() { + std::unordered_map rsum_map; + const Opm::SummaryConfig summaryConfig = simulator_.vanguard().summaryConfig(); + + rsum_map.emplace("FIPNUM", makeRegionSum(this->regions_.at("FIPNUM"), true)); + for (const auto& fip_region : summaryConfig.fip_regions()) { + if (fip_region == "FIPNUM") + continue; + + const auto& region = this->regions_.at(fip_region); + rsum_map.emplace(fip_region, makeRegionSum(region, false)); + } + return rsum_map; + } + + Scalar sum(const ScalarBuffer& v) { + return std::accumulate(v.begin(), v.end(), Scalar{0}); + } + + void updateSummaryRegionValues(const std::unordered_map& rsum_map, + std::map& miscSummaryData, + std::map>& regionData) const { + + const Opm::SummaryConfig summaryConfig = simulator_.vanguard().summaryConfig(); + + // The field summary vectors should only use the FIPNUM based region sum. + { + const auto& rsum = rsum_map.at("FIPNUM"); + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) { + std::string key = "F" + FipDataType::EclString(fip_type); + if (summaryConfig.hasKeyword(key)) + miscSummaryData[key] = rsum.fieldFipValues[fip_type]; + } + if (summaryConfig.hasKeyword("FOE") && this->fieldInitialInplace_) + miscSummaryData["FOE"] = rsum.fieldFipValues[FipDataType::OilInPlace] + / this->fieldInitialInplace_.value()[FipDataType::OilInPlace]; + + if (summaryConfig.hasKeyword("FPR")) + miscSummaryData["FPR"] = pressureAverage_(rsum.fieldPressurePvHydrocarbon, + rsum.fieldPvHydrocarbon, + rsum.fieldPressurePv, + rsum.fieldFipValues[FipDataType::PoreVolume], + true); + + if (summaryConfig.hasKeyword("FPRP")) + miscSummaryData["FPRP"] = pressureAverage_(rsum.fieldPressurePvHydrocarbon, + rsum.fieldPvHydrocarbon, + rsum.fieldPressurePv, + rsum.fieldFipValues[FipDataType::PoreVolume], + false); + + } + + // The region summary vectors should loop through the FIPxxx regions to + // support the RPR__xxx summary keywords. + { + for (int fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) { + for (const auto& node : this->regionNodes_[fip_type]) { + const auto& rsum = rsum_map.at(node.fip_region()); + regionData[node.keyword()] = rsum.regFipValues[fip_type]; + } + } + + // The exact same quantity is calculated for RPR and RPRP - is that correct? + for (const auto& node : this->RPRNodes_) { + const auto& rsum = rsum_map.at(node.fip_region()); + regionData[node.keyword()] = pressureAverage_(rsum.regPressurePvHydrocarbon, + rsum.regPvHydrocarbon, + rsum.regPressurePv, + rsum.regFipValues[FipDataType::PoreVolume], + true); + } + + for (const auto& node : this->RPRPNodes_) { + const auto& rsum = rsum_map.at(node.fip_region()); + regionData[node.keyword()] = pressureAverage_(rsum.regPressurePvHydrocarbon, + rsum.regPvHydrocarbon, + rsum.regPressurePv, + rsum.regFipValues[FipDataType::PoreVolume], + false); + } + } + } + + + void outputFipLogImpl(const RegionSum& rsum) const { + + + { + Scalar fieldHydroCarbonPoreVolumeAveragedPressure = pressureAverage_(rsum.fieldPressurePvHydrocarbon, + rsum.fieldPvHydrocarbon, + rsum.fieldPressurePv, + rsum.fieldFipValues[FipDataType::PoreVolume], + true); + std::array initial_values = *this->fieldInitialInplace_; + std::array current_values = rsum.fieldFipValues; + fipUnitConvert_(initial_values); + fipUnitConvert_(current_values); + + pressureUnitConvert_(fieldHydroCarbonPoreVolumeAveragedPressure); + outputRegionFluidInPlace_(initial_values, + current_values, + fieldHydroCarbonPoreVolumeAveragedPressure); + } + + for (size_t reg = 0; reg < rsum.ntFip; ++reg) { + std::array initial_values; + std::array current_values; + + for (std::size_t fip_type = 0; fip_type < FipDataType::numFipTypes; fip_type++) { + initial_values[fip_type] = this->regionInitialInplace_.value()[fip_type][reg]; + current_values[fip_type] = rsum.regFipValues[fip_type][reg]; + } + fipUnitConvert_(initial_values); + fipUnitConvert_(current_values); + + Scalar regHydroCarbonPoreVolumeAveragedPressure + = pressureAverage_(rsum.regPressurePvHydrocarbon[reg], + rsum.regPvHydrocarbon[reg], + rsum.regPressurePv[reg], + rsum.regFipValues[FipDataType::PoreVolume][reg], + true); + pressureUnitConvert_(regHydroCarbonPoreVolumeAveragedPressure); + outputRegionFluidInPlace_(initial_values, current_values, regHydroCarbonPoreVolumeAveragedPressure, reg + 1); + } + } + + + + // write Fluid In Place to output log void outputFipLog(std::map& miscSummaryData, std::map>& regionData, const bool substep) { - const auto& comm = simulator_.gridView().comm(); - auto maxElement = std::max_element(fipnum_.begin(), fipnum_.end()); - size_t ntFip = 0; - if ( maxElement != fipnum_.end() ) { - ntFip = *maxElement; - } - ntFip = comm.max(ntFip); + auto rsum_map = this->accumulateRegionSums(); + if (!isIORank_()) + return; - // sum values over each region - ScalarBuffer regionFipValues[FipDataType::numFipValues]; - for (int i = 0; i < FipDataType::numFipValues; i++) { - regionFipValues[i] = computeFipForRegions_(fip_[i], fipnum_, ntFip); - if (isIORank_() && origRegionValues_[i].empty()) - origRegionValues_[i] = regionFipValues[i]; - } - - // sum all region values to compute the field total - std::vector fieldNum(ntFip, 1); - ScalarBuffer fieldFipValues(FipDataType::numFipValues, 0.0); - bool comunicateSum = false; // the regionValues are already summed over all ranks. - for (int i = 0; i& region) { - const auto& gridView = simulator_.vanguard().gridView(); - if (simulator_.vanguard().eclState().fieldProps().has_int("FIPNUM")) - fipnum_ = simulator_.vanguard().eclState().fieldProps().get_int("FIPNUM"); - else - fipnum_.resize(gridView.size(/*codim=*/0), 1); - ElementContext elemCtx(simulator_); ElementIterator elemIt = simulator_.gridView().template begin(); const ElementIterator& elemEndIt = simulator_.gridView().template end(); @@ -1787,19 +1926,19 @@ private: for (; elemIt != elemEndIt; ++elemIt, ++elemIdx) { const Element& elem = *elemIt; if (elem.partitionType() != Dune::InteriorEntity) - fipnum_[elemIdx] = 0; + region[elemIdx] = 0; } } // Sum Fip values over regions. - ScalarBuffer computeFipForRegions_(const ScalarBuffer& fip, std::vector& regionId, size_t maxNumberOfRegions, bool commSum = true) + ScalarBuffer regionSum(const ScalarBuffer& property, const std::vector& regionId, size_t maxNumberOfRegions) const { ScalarBuffer totals(maxNumberOfRegions, 0.0); - if (fip.empty()) + if (property.empty()) return totals; - assert(regionId.size() == fip.size()); + assert(regionId.size() == property.size()); for (size_t j = 0; j < regionId.size(); ++j) { const int regionIdx = regionId[j] - 1; // the cell is not attributed to any region. ignore it! @@ -1807,18 +1946,17 @@ private: continue; assert(regionIdx < static_cast(maxNumberOfRegions)); - totals[regionIdx] += fip[j]; - } - if (commSum) { - const auto& comm = simulator_.gridView().comm(); - for (size_t i = 0; i < maxNumberOfRegions; ++i) - totals[i] = comm.sum(totals[i]); + totals[regionIdx] += property[j]; } + const auto& comm = simulator_.gridView().comm(); + for (size_t i = 0; i < maxNumberOfRegions; ++i) + totals[i] = comm.sum(totals[i]); + return totals; } - ScalarBuffer pressureAverage_(const ScalarBuffer& pressurePvHydrocarbon, const ScalarBuffer& pvHydrocarbon, const ScalarBuffer& pressurePv, const ScalarBuffer& pv, bool hydrocarbon) + ScalarBuffer pressureAverage_(const ScalarBuffer& pressurePvHydrocarbon, const ScalarBuffer& pvHydrocarbon, const ScalarBuffer& pressurePv, const ScalarBuffer& pv, bool hydrocarbon) const { size_t size = pressurePvHydrocarbon.size(); assert(pvHydrocarbon.size() == size); @@ -1832,7 +1970,7 @@ private: return fraction; } - Scalar pressureAverage_(const Scalar& pressurePvHydrocarbon, const Scalar& pvHydrocarbon, const Scalar& pressurePv, const Scalar& pv, bool hydrocarbon) + Scalar pressureAverage_(const Scalar& pressurePvHydrocarbon, const Scalar& pvHydrocarbon, const Scalar& pressurePv, const Scalar& pv, bool hydrocarbon) const { if (pvHydrocarbon > 1e-10 && hydrocarbon) return pressurePvHydrocarbon / pvHydrocarbon; @@ -1840,7 +1978,7 @@ private: return pressurePv / pv; } - void fipUnitConvert_(ScalarBuffer& fip) + void fipUnitConvert_(std::array& fip) const { const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { @@ -1872,7 +2010,7 @@ private: } } - void pressureUnitConvert_(Scalar& pav) + void pressureUnitConvert_(Scalar& pav) const { const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { @@ -1890,7 +2028,9 @@ private: } } - void outputRegionFluidInPlace_(const ScalarBuffer& oip, const ScalarBuffer& cip, const Scalar& pav, const int reg) + void outputRegionFluidInPlace_(const std::array& oip, + const std::array& cip, + const Scalar& pav, const int reg = 0) const { if (forceDisableFipOutput_) return; @@ -1901,7 +2041,7 @@ private: const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); std::ostringstream ss; - if (!reg) { + if (reg == 0) { ss << " ===================================================\n" << " : Field Totals :\n"; } @@ -2043,22 +2183,6 @@ private: Opm::OpmLog::note(ss.str()); } - std::string fipEnumToString_(int i) - { - typedef typename FipDataType::FipId FipId; - switch(static_cast(i)) - { - case FipDataType::WaterInPlace: return "WIP"; - case FipDataType::OilInPlace: return "OIP"; - case FipDataType::GasInPlace: return "GIP"; - case FipDataType::OilInPlaceInLiquidPhase: return "OIPL"; - case FipDataType::OilInPlaceInGasPhase: return "OIPG"; - case FipDataType::GasInPlaceInLiquidPhase: return "GIPL"; - case FipDataType::GasInPlaceInGasPhase: return "GIPG"; - case FipDataType::PoreVolume: return "PV"; - } - return "ERROR"; - } std::string WPEnumToString_(int i) { typedef typename WellProdDataType::WPId WPId; @@ -2171,10 +2295,17 @@ private: std::vector failedCellsPb_; std::vector failedCellsPd_; - std::vector fipnum_; - ScalarBuffer fip_[FipDataType::numFipValues]; - ScalarBuffer origTotalValues_; - ScalarBuffer origRegionValues_[FipDataType::numFipValues]; + std::unordered_map> regions_; + std::array fip_; + std::optional> regionInitialInplace_; + std::optional> fieldInitialInplace_; + + std::vector RPRNodes_; + std::vector RPRPNodes_; + std::array, FipDataType::numFipTypes> regionNodes_; + + + ScalarBuffer hydrocarbonPoreVolume_; ScalarBuffer pressureTimesPoreVolume_; ScalarBuffer pressureTimesHydrocarbonVolume_;