Merge pull request #3397 from bska/report-dynpv

Report Pressure Dependent Pore Volume in PRT File

ebos/eclgenericoutputblackoilmodule.cc

@@ -38,6 +38,7 @@
 #include <opm/parser/eclipse/Units/Units.hpp>
 
 #include <cassert>
+#include <initializer_list>
 #include <iomanip>
 #include <sstream>
 #include <stdexcept>
@@ -827,33 +828,41 @@ doAllocBuffers(unsigned bufferSize,
         }
     }
 
-    outputFipRestart_ = false;
+    this->outputFipRestart_ = false;
-    computeFip_ = 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;
-        }
-        else
-            fip_[phase].clear();
             }
 
-    if (!substep || summaryConfig_.hasKeyword("FPR") || summaryConfig_.hasKeyword("FPRP") || !this->RPRNodes_.empty()) {
+            this->fip_[phase].resize(bufferSize, 0.0);
-        fip_[Inplace::Phase::PoreVolume].resize(bufferSize, 0.0);
+            this->computeFip_ = true;
-        hydrocarbonPoreVolume_.resize(bufferSize, 0.0);
-        pressureTimesPoreVolume_.resize(bufferSize, 0.0);
-        pressureTimesHydrocarbonVolume_.resize(bufferSize, 0.0);
         }
         else {
-        hydrocarbonPoreVolume_.clear();
+            this->fip_[phase].clear();
-        pressureTimesPoreVolume_.clear();
+        }
-        pressureTimesHydrocarbonVolume_.clear();
+    }
+
+    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 {
+        this->dynamicPoreVolume_.clear();
+        this->hydrocarbonPoreVolume_.clear();
+        this->pressureTimesPoreVolume_.clear();
+        this->pressureTimesHydrocarbonVolume_.clear();
     }
 
     // Well RFT data
@@ -1074,32 +1083,25 @@ void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
 fipUnitConvert_(std::unordered_map<Inplace::Phase, Scalar>& fip) const
 {
     const UnitSystem& units = eclState_.getUnits();
-    if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
+    using M = UnitSystem::measure;
-        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);
+    const auto unit_map = std::unordered_map<Inplace::Phase, M> {
-        fip[Inplace::Phase::OilInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::OilInLiquidPhase], unit::stb);
+        {Inplace::Phase::WATER,             M::liquid_surface_volume},
-        fip[Inplace::Phase::OilInGasPhase] = unit::convert::to(fip[Inplace::Phase::OilInGasPhase], unit::stb);
+        {Inplace::Phase::OIL,               M::liquid_surface_volume},
-        fip[Inplace::Phase::GAS] = unit::convert::to(fip[Inplace::Phase::GAS], 1000*unit::cubic(unit::feet));
+        {Inplace::Phase::OilInLiquidPhase,  M::liquid_surface_volume},
-        fip[Inplace::Phase::GasInLiquidPhase] = unit::convert::to(fip[Inplace::Phase::GasInLiquidPhase], 1000*unit::cubic(unit::feet));
+        {Inplace::Phase::OilInGasPhase,     M::liquid_surface_volume},
-        fip[Inplace::Phase::GasInGasPhase] = unit::convert::to(fip[Inplace::Phase::GasInGasPhase], 1000*unit::cubic(unit::feet));
+        {Inplace::Phase::GAS,               M::gas_surface_volume},
-        fip[Inplace::Phase::PoreVolume] = unit::convert::to(fip[Inplace::Phase::PoreVolume], unit::stb);
+        {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);
         }
-    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.");
     }
 }
 
@@ -1107,20 +1109,8 @@ template<class FluidSystem, class Scalar>
 void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
 pressureUnitConvert_(Scalar& pav) const
 {
-    const UnitSystem& units = eclState_.getUnits();
+    pav = this->eclState_.getUnits()
-    if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
+        .from_si(UnitSystem::measure::pressure, pav);
-        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.");
-    }
 }
 
 template<class FluidSystem, class Scalar>
@@ -1133,11 +1123,25 @@ outputRegionFluidInPlace_(std::unordered_map<Inplace::Phase, Scalar> 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";
@@ -1307,9 +1311,9 @@ pressureAverage_(const Scalar& pressurePvHydrocarbon,
                  const Scalar& pvHydrocarbon,
                  const Scalar& pressurePv,
                  const Scalar& pv,
-             bool hydrocarbon)
+                 const bool    hydrocarbon)
 {
-    if (pvHydrocarbon > 1e-10 && hydrocarbon)
+    if (hydrocarbon && (pvHydrocarbon > 1e-10))
         return pressurePvHydrocarbon / pvHydrocarbon;
 
     return pressurePv / pv;
@@ -1322,17 +1326,22 @@ pressureAverage_(const ScalarBuffer& pressurePvHydrocarbon,
                  const ScalarBuffer& pvHydrocarbon,
                  const ScalarBuffer& pressurePv,
                  const ScalarBuffer& pv,
-             bool hydrocarbon)
+                 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) {
+    for (std::size_t i = 0; i < size; ++i) {
-        fraction[i] = pressureAverage_(pressurePvHydrocarbon[i], pvHydrocarbon[i], pressurePv[i], pv[i], hydrocarbon);
+        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,
+        outputRegionFluidInPlace_(std::move(initial_values),
-                                  current_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<class FluidSystem,class Scalar>
 void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
 update(Inplace& inplace,
        const std::string& region_name,
-       Inplace::Phase phase,
+       const Inplace::Phase phase,
-       std::size_t ntFip,
+       const std::size_t ntFip,
-       const std::vector<double>& values)
+       const ScalarBuffer& values)
 {
-    double sum = 0;
+    double sum = 0.0;
-    for (std::size_t region_number = 0; region_number < ntFip; region_number++) {
+    for (std::size_t region_number = 0; region_number < ntFip; ++region_number) {
-        inplace.add( region_name, phase, region_number + 1, values[region_number] );
+        const auto rval = static_cast<double>(values[region_number]);
-        sum += values[region_number];
+        inplace.add(region_name, phase, region_number + 1, rval);
+        sum += rval;
     }
-    inplace.add( phase, sum );
+    inplace.add(phase, sum);
 }
 
 template<class FluidSystem,class Scalar>
 void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
 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));
+    auto update_inplace =
-    update(inplace, region_name, Inplace::Phase::HydroCarbonPV, ntFip, this->regionSum(this->hydrocarbonPoreVolume_, region, ntFip, comm));
+        [&inplace, &region, &region_name, &comm, ntFip, this]
-    update(inplace, region_name, Inplace::Phase::PressureHydroCarbonPV, ntFip, this->regionSum(this->pressureTimesHydrocarbonVolume_, region, ntFip, comm));
+        (const Inplace::Phase       phase,
+         const std::vector<Scalar>& value)
+    {
+        update(inplace, region_name, phase, ntFip,
+               this->regionSum(value, region, ntFip, comm));
+    };
 
-    for (const auto& phase : Inplace::phases())
+    update_inplace(Inplace::Phase::PressurePV,
-        update(inplace, region_name, phase, ntFip, this->regionSum(this->fip_[phase], region, ntFip, comm));
+                   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<class FluidSystem,class Scalar>
@@ -1501,9 +1540,8 @@ accumulateRegionSums(const Comm& comm)
 {
     Inplace inplace;
 
-    for (const auto& [region_name, _] : this->regions_) {
+    for (const auto& region : this->regions_) {
-        (void)_;
+        makeRegionSum(inplace, region.first, comm);
-        makeRegionSum(inplace, region_name, comm);
     }
 
     // The first time the outputFipLog function is run we store the inplace values in
@@ -1537,55 +1575,69 @@ 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);
+            const std::string key = "F" + EclString(phase);
-            if (summaryConfig_.hasKeyword(key))
+            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"))
+        if (this->summaryConfig_.hasKeyword("FPR")) {
-            miscSummaryData["FPR"] = pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV),
+            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("FPRP")) {
-        if (summaryConfig_.hasKeyword("FPRP"))
+            miscSummaryData["FPRP"] =
-            miscSummaryData["FPRP"] = pressureAverage_(inplace.get(Inplace::Phase::PressureHydroCarbonPV),
+                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_) {
-        for (const auto& node : this->RPRNodes_)
+            regionData[node.keyword()] =
-            regionData[node.keyword()] = pressureAverage_(inplace.get_vector(node.fip_region(), Inplace::Phase::PressureHydroCarbonPV),
+            pressureAverage_(get_vector(node, Inplace::Phase::PressureHydroCarbonPV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::HydroCarbonPV),
+                             get_vector(node, Inplace::Phase::HydroCarbonPV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::PressurePV),
+                             get_vector(node, Inplace::Phase::PressurePV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::PoreVolume),
+                             get_vector(node, Inplace::Phase::PoreVolume),
                              true);
+        }
 
-
+        for (const auto& node : this->RPRPNodes_) {
-        for (const auto& node : this->RPRPNodes_)
+            regionData[node.keyword()] =
-            regionData[node.keyword()] = pressureAverage_(inplace.get_vector(node.fip_region(), Inplace::Phase::PressureHydroCarbonPV),
+            pressureAverage_(get_vector(node, Inplace::Phase::PressureHydroCarbonPV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::HydroCarbonPV),
+                             get_vector(node, Inplace::Phase::HydroCarbonPV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::PressurePV),
+                             get_vector(node, Inplace::Phase::PressurePV),
-                                                          inplace.get_vector(node.fip_region(), Inplace::Phase::PoreVolume),
+                             get_vector(node, Inplace::Phase::PoreVolume),
                              false);
         }
+    }
 }
 
 template class EclGenericOutputBlackoilModule<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,double>;

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<int>& regionId,
-                                  size_t maxNumberOfRegions,
+                                  const std::size_t maxNumberOfRegions,
                                   const Comm& comm);
 
     static int regionMax(const std::vector<int>& region,
@@ -293,9 +293,9 @@ protected:
 
     static void update(Inplace& inplace,
                        const std::string& region_name,
-                       Inplace::Phase phase,
+                       const Inplace::Phase phase,
-                       std::size_t ntFip,
+                       const std::size_t ntFip,
-                       const std::vector<double>& values);
+                       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_;

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 =
+        const auto totVolume =
-            elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
+            elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
-            * intQuants.porosity().value();
+        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))

opm/simulators/wells/BlackoilWellModel.hpp

@@ -222,7 +222,8 @@ namespace Opm {
 
             data::Wells wellData() const
             {
-                auto wsrpt = this->wellState().report(UgGridHelpers::globalCell(grid()),
+                auto wsrpt = this->wellState()
+                    .report(UgGridHelpers::globalCell(grid()),
                             [this](const int well_ndex) -> bool
                 {
                     return this->wasDynamicallyShutThisTimeStep(well_ndex);